Jako vstupní zdroje dat byla použita již transformovaná data od zdravotních pojišťoven ČPZP a OZP upravená do datábázového formátu a uložena v lokální DuckDB databázi. Vstupní soubory byly očištěny o předpisy skupiny L04, které by eventuelně mohly sloužit k samostatné analýze, ale nedají se jednoznačně označit jako substituty kortikoidsteroidů.

Během vyhodnocování byli vždy vybráni pouze ti klienti, u kterých máme jistotu, že jsou pojištěni po celé sledované období. Pokud by noví klienti nebyli vyfiltrováni, hrozilo by výrazné zkreslení prvopředpisů a zkreslení celkového trendu způsobené fluktuací. Pokud klient v průběhu sledování zemřel, je bráno datum úmrtí jako datum konce pojištění. V případě dlouhodobého vyhodnocování jsou tedy vybráni klienti s počátkem pojištění před rokem 2015 a bez konce pojištění. Pro krátkodobé vyhodnocení (SCCS) byli vybráni klienti, kteří měli platné pojištění minimálně +- 365 dní od rozhodného data (vrchol vakcinace dané kohorty plus jeden týden).

Věk klienta je počítán k roku 2021, pokud není uvedeno jinak (např. U výpočtu rizika prvopředpisu podle věku).

Jako očkovaní klienti byli vybráni ti s alespoň jedním záznamem o vakcinaci a jako datum vakcinace bylo zvoleno datum první dávky.

Předpisy injekcí / infuzí byly agregovány do jednoho předpisu pouze u krátkodobé analýzy (SCCS). V ostatních případech sledujeme buď celkový trend (kde je výkyv v počtu předpisů způsobený větším množstvím injekcí také jeden z možných výsledků) nebo prvopředpisy (kde je množství irelevantní).

Při posuzování časových řad pomocí ETS modelu (příp. pomocí STL s lineární regresí) je posuzován nastavený trend do roku 2020 proti reálnému průběhu po tomto roce. Cíl tohoto přístupu je identifikovat, jestli měl COVID-19 nějaký vliv na celkový trend a případně, v jaké skupině dochází k největší odchylce.

V rámci analýz dlouhodobého dopadu jsou porovnávána data 2018-2020 proti datům z let 2022 a 2023. Smysl tohoto přístupu je odhalit, jestli měl COVID-19 (a potažmo očkování) nějaký vliv na období po relativním uklidnění. Data z let 2020 a 2021 jsou také vyřazena z důvodu nestability a výrazných sezónních propadů.

Pro dlouhodobé dopady byly použity 3 metriky, které je možné i částečně vyhodnotit vedle sebe. Vždy je bráno v úvahu, že pokud dochází k měření určitého stavu, musí mít všichni stejnou příležitost k dosažení daného stavu. Pro konkrétní metriky tedy:

  • Prvopředpisovost: Data z počátku měření vykazují jako prvopředpis v zásadě každý předpis. To je samozřejmě chybný odhad a tedy jsou ze všech prvopředpisových metrik vyjmuty roky 2015 a 2016. Zároveň pokud je míra prvopředpisovosti stále identická, tak celkový úhrn prvopředpisů v čase klesá (nikdo nemůže dostat prvopředpis dvakrát). Z toho důvodu je prvopředpisovost vždy počítaná jako počet prvopředpisů na 1000 dosud nepředepisovaných osob.

  • Multipředpisovost: Namísto celkového úhrnu předpisů je počítán poměr osob s více než jedním předpisem za zvolené období (měsíc nebo rok) na 1000 předepisovaných osob. Tento přístup umožní nejen vysledovat samotnou míru zvýšení či snížení, ale i odhadnout o jaký typ léčby jde.

Smysl by dávalo měřit i celkový počet předpisů na 1000 osob všeobecně a počty unikátních osob v daném období. Na obojí má ale výrazný vliv retence (tedy stav, kdy předepsaný klient dostává předpisy i ve všech dalších sledovaných časových intervalech) a vzhledem k tomu, že analýza je postavená na uzavřené skupině, kde nedochází k žádné fluktuaci klientů, tak by bylo nutné obě tyto metriky očistit o přirozenou míru udržení (kterou nelze jednoznačně odhadnout) a celkově by takové měření bylo náročné na interpretaci.

Krátkodobá analýza aplikuje metodu SCCS na prvopředpisy a veškeré předpisy. Smysl tohoto přístupu je odhalit bezprostřední rizika vakcinace a z pohledu celé analýzy je tento postup nejdůležitějším pro posouzení i z toho důvodu, že v kratším časovém horizontu lze předpokládat nižší vliv zavádějících proměnných způsobených časem.

Kontrolní součty z pohledu dlouhodobého sledování:

  • Počet očkovaných jedinců u jednotlivých pojišťoven: OZP 433 923 ČPZP 646 862
  • Počet neočkovaných jedinců bez předpisu u jednotlivých pojišťoven: OZP 106 160 ČPZP 248 729
  • Počet předpisů kortikoidů ve věkové skupině 31-50 let: 671 453
  • Počet očkovaných 31-50 let: 358 346
  • Počet předpisů pro očkované 31-50 let pro roky 2020-2022: 251 194 (včetně L04), 176 973 (H02)
  • Počet neplatných záznamů: Neplatné záznamy jsou filtrovány samostaně na úrovni výpočtů a SQL dotazů. Například při analýze formy podání bylo nalezeno 56 záznamů s datem předpisu, ale s prázdnou formou podání.
  • Počet prvopředpisů pro roky 2020-2022: 145 195 (včetně L04), 144 793 (H02)

1 Celkový přehled

1.1 Předpisy

prescriptions <- dbGetQuery(con,
                          "
                          SELECT
                            c.client_id,
                            c.n_vacc,
                            c.n_presc,
                            c.sex,
                            c.birthyear,
                            2021 - c.birthyear AS age_2021,
                            FIRST(p.date) OVER (PARTITION BY c.client_id ORDER BY p.presc_order ASC) AS first_presc_date,
                            p.date AS presc_date,
                            v.date AS vacc_date,
                            d.drug_form,
                            d.ATC_group,
                            d.Prednison_equiv AS prednison_equiv,
                            d.force AS strength,
                            'CPZP' AS ins_company
                          FROM
                            cpzp_clients c
                          LEFT JOIN cpzp_ins_start_end i ON c.client_id = i.client_id
                          LEFT JOIN cpzp_prescriptions p ON c.client_id = p.client_id
                          LEFT JOIN cpzp_vaccinations v ON c.client_id = v.client_id AND event_order = 1
                          LEFT JOIN cpzp_drugs_catalog d ON p.presc_id  = d.presc_id 
                          WHERE
                            i.ins_start_date < '2015-01-01' AND i.ins_end_date = '9999-12-31'

                            
                          UNION
                          
                          SELECT
                            c.client_id,
                            c.n_vacc,
                            c.n_presc,
                            c.sex,
                            c.birthyear,
                            2021 - c.birthyear AS age_2021,
                            FIRST(p.date) OVER (PARTITION BY c.client_id ORDER BY p.date ASC) AS first_presc_date,
                            p.date AS presc_date,
                            v.date AS vacc_date,
                            d.drug_form,
                            d.ATC_group,
                            d.Prednison_equiv AS prednison_equiv,
                            d.strength,
                            'OZP' AS ins_company
                          FROM
                            ozp_clients c
                          LEFT JOIN ozp_ins_start_end i ON c.client_id = i.client_id
                          LEFT JOIN ozp_prescriptions p ON c.client_id = p.client_id
                          LEFT JOIN ozp_vaccinations v ON c.client_id = v.client_id AND event_order = 1
                          LEFT JOIN ozp_drugs_catalog d ON p.presc_id = d.presc_id 
                          WHERE
                            i.ins_start_date < '2015-01-01' AND i.ins_end_date IS NULL
                          ") %>%
  mutate(client_id = as.numeric(factor(client_id)),
         vacc_yearmonth = floor_date(vacc_date, unit = "month"),
         presc_yearmonth = floor_date(presc_date, unit = "month"),
         vaccinated = ifelse(n_vacc > 0, "Vakcinace", "Bez vakcinace"),
         sex = ifelse(sex == "Z", "F", sex),
         first_presc_flag = ifelse(!is.na(presc_date) & first_presc_date == presc_date,1,0)) %>%
mutate(age_class = case_when(
    age_2021 >= 5 & age_2021 <= 11 ~ "5-11",
    age_2021 >= 12 & age_2021 <= 15 ~ "12-15",
    age_2021 >= 16 & age_2021 <= 29 ~ "16-29",
    age_2021 >= 30 & age_2021 <= 34 ~ "30-34",
    age_2021 >= 35 & age_2021 <= 39 ~ "35-39",
    age_2021 >= 40 & age_2021 <= 44 ~ "40-44",
    age_2021 >= 45 & age_2021 <= 49 ~ "45-49",
    age_2021 >= 50 & age_2021 <= 54 ~ "50-54",
    age_2021 >= 55 & age_2021 <= 59 ~ "55-59",
    age_2021 >= 60 & age_2021 <= 64 ~ "60-64",
    age_2021 >= 65 & age_2021 <= 69 ~ "65-69",
    age_2021 >= 70 & age_2021 <= 79 ~ "70-79",
    age_2021 >= 80 ~ "80+",
    TRUE ~ NA_character_
  ),
  age_class = factor(age_class, levels = c(
    "5-11", "12-15", "16-29", "30-34", "35-39", "40-44",
    "45-49", "50-54", "55-59", "60-64", "65-69", "70-79", "80+"
  ))
)

1.1.1 Vývoj počtu předpisů v čase

Na první pohled je viditelná odchylka od trendu během roku 2020 a následný strmější nárůst. Zároveň je také zřejmá vyšší míra předpisovosti u žen. Nižší počet předpisů u neočkované populace je jednoduše vysvětlitelný nižším počtem neočkovaných a pravděpodobně i jiným přístupem k prevenci.

total_presc_by_month <- prescriptions %>%
  filter(!is.na(presc_yearmonth)) %>%
  group_by(presc_yearmonth) %>%
  summarise(total_prescriptions = n()) %>%
  ungroup()

ggplot(total_presc_by_month, aes(x = presc_yearmonth, y = total_prescriptions)) +
  geom_line(color = "steelblue", linewidth = 1, alpha = 0.6) +
  stat_smooth(method = "loess", span = 0.3, se = TRUE, linewidth = 1.5, color = "firebrick") +
  scale_x_date(
    date_labels = "%Y",
    date_breaks = "1 year",
    expand = expansion(mult = c(0.01, 0.01))
  ) +
  labs(
    title = "Celkový úhrn předpisů",
    x = "Rok",
    y = "Počet předpisů"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    axis.text.x = element_text(color = "black", size = 12, face = "bold"),
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5)
  )

grouped_presc <- prescriptions %>%
  filter(!is.na(presc_yearmonth)) %>%
  group_by(presc_yearmonth, sex, vaccinated) %>%
  summarise(total_prescriptions = n()) %>%
  ungroup()

ggplot(grouped_presc, aes(x = presc_yearmonth, y = total_prescriptions,
                         color = vaccinated, linetype = sex)) +
  geom_vline(xintercept = as.numeric(as.Date("2020-12-01")), linetype = "dotted", color = "gray70", linewidth = 0.4) +
  geom_line(linewidth = 1, alpha = 0.8) +
  stat_smooth(method = "loess", span = 0.3, se = TRUE, linewidth = 1.5, aes(color = vaccinated), linetype = "solid") +
  scale_x_date(
    date_labels = "%Y",
    date_breaks = "1 year",
    expand = expansion(mult = c(0.01, 0.01))
  ) +
  scale_color_manual(values = c("Vakcinace" = "darkgreen", "Bez vakcinace" = "firebrick")) +
  labs(
    title = "Předpisy podle pohlaví a vakcinace",
    x = "Rok",
    y = "Počet předpisů",
    color = "Vakcinace",
    linetype = "Pohlaví"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    axis.text.x = element_text(color = "black", size = 12, face = "bold"),
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5)
  )

1.1.2 Riziko prvopředpisu podle věku

Zde ze zřejmých důvodů není použit věk v roce 2021, ale věk v roce prvopředpisu. Riziko je počítané jako podíl prvopředpisů na počet osob v daném věku, které ještě žádný předpis neměly. Vyřazené jsou prvopředpisy v roce 2015 a 2016, které budou pravděpodobně vykazovat vysokou chybovost u lidí předepisovaných sezónně. Distribuce z logických důvodů (nemáme k dispozici veškerá historická data pro každého klienta) částečně kopíruje populační křivku, výsledek je tedy nutné brát s jistou rezervou.

first_presc_dynamic_age <- dbGetQuery(con,
                          "
                          SELECT
                            c.client_id,
                            c.n_vacc,
                            c.n_presc,
                            c.sex,
                            c.birthyear,
                            FIRST(p.date) OVER (PARTITION BY c.client_id ORDER BY p.date ASC) AS first_presc_date,
                            'CPZP' AS ins_company
                          FROM
                            cpzp_clients c
                          LEFT JOIN cpzp_ins_start_end i ON c.client_id = i.client_id
                          LEFT JOIN cpzp_prescriptions p ON c.client_id = p.client_id AND p.presc_order = 1
                          WHERE
                            i.ins_start_date < '2015-01-01' AND i.ins_end_date = '9999-12-31'
                            
                          UNION
                          
                          SELECT
                            c.client_id,
                            c.n_vacc,
                            c.n_presc,
                            c.sex,
                            c.birthyear,
                            FIRST(p.date) OVER (PARTITION BY c.client_id ORDER BY p.date ASC) AS first_presc_date,
                            'OZP' AS ins_company
                          FROM
                            ozp_clients c
                          LEFT JOIN ozp_ins_start_end i ON c.client_id = i.client_id
                          LEFT JOIN ozp_prescriptions p ON c.client_id = p.client_id AND p.presc_order = 1
                          WHERE
                            i.ins_start_date < '2015-01-01' AND i.ins_end_date IS NULL
                          ") %>%
  mutate(first_presc_yearmonth = floor_date(first_presc_date, unit = "month"),
         sex = ifelse(sex == "Z", "F", sex),
         vacc_status = ifelse(n_vacc > 0, "Vakcinace", "Bez vakcinace"),
         age_at_presc = year(first_presc_date) - birthyear)

obs_start <- 2017
obs_end <- 2023
max_age <- 80


filtered_data <- first_presc_dynamic_age %>%
  mutate(
    year_of_presc = year(first_presc_date),
    age_of_presc = year_of_presc - birthyear,
    age_at_start = obs_start - birthyear,
    age_at_end = obs_end - birthyear
  ) %>%
  filter(age_at_end >= 0) %>%   
  filter(age_at_start <= max_age)

age_seq <- 0:max_age

at_risk_by_age <- lapply(age_seq, function(a) {
  n <- filtered_data %>%
    filter(is.na(age_of_presc) | age_of_presc > a) %>%
    n_distinct(.$client_id)
  data.frame(age = a, n_at_risk = n)
}) %>%
  bind_rows()

first_presc_by_age <- filtered_data %>%
  filter(!is.na(age_of_presc), year(first_presc_date) >= obs_start) %>%
  count(age = age_of_presc, name = "n_first_presc") %>%
  filter(age <= max_age)

age_analysis <- left_join(at_risk_by_age, first_presc_by_age, by = "age") %>%
  mutate(
    n_first_presc = replace_na(n_first_presc, 0),
    prop_first_presc = n_first_presc / n_at_risk
  )

ggplot(age_analysis, aes(x = age, y = prop_first_presc * 1000)) +
  geom_line(color = "steelblue") +
  geom_point(aes(color = n_first_presc)) +
  labs(
    title = "Riziko prvopředpisu podle věku",
    x = "Věk",
    y = "Počet prvopředpisů na 1000 osob bez předpisu",
    color = "Počet prvopředpisů"
  ) +
  theme_minimal()

1.2 Vakcinace

assign_age_class <- function(birthyear) {
  age <- 2021 - birthyear
  
  cut(
    age,
    
    breaks = c(-Inf, 4, 11, 15, 29, 34, 39, 44, 49, 54, 59, 64, 69, 79, Inf),
    labels = c(
      "<5",
      "5-11",
      "12-15",
      "16-29",
      "30-34",
      "35-39",
      "40-44",
      "45-49",
      "50-54",
      "55-59",
      "60-64",
      "65-69",
      "70-79",
      "80+"
    ),
    right = TRUE,
    ordered_result = TRUE
  )
}

vaccination_first_shot_by_birthyear <- dbGetQuery(con,
                                "
                                SELECT DISTINCT
                                v.date,
                                c.birthyear,
                                COUNT(DISTINCT v.client_id) OVER (
                                PARTITION BY v.date, c.birthyear) AS vaccinated_total_in_day,
                                'OZP' AS ins_company
                                FROM ozp_vaccinations v
                                LEFT JOIN ozp_clients c ON v.client_id = c.client_id
                                LEFT JOIN ozp_ins_start_end i ON c.client_id = i.client_id
                                WHERE
                                v.event_order = 1
                                AND i.ins_start_date < '2020-01-01' AND i.ins_end_date > '2023-01-01'
                                
                                UNION
                                
                                SELECT DISTINCT
                                v.date,
                                c.birthyear,
                                COUNT(DISTINCT v.client_id) OVER (
                                PARTITION BY v.date, c.birthyear) AS vaccinated_total_in_day,
                                'CPZP' AS ins_company
                                FROM cpzp_vaccinations v
                                LEFT JOIN cpzp_clients c ON v.client_id = c.client_id
                                LEFT JOIN cpzp_ins_start_end i ON c.client_id = i.client_id
                                WHERE
                                v.event_order = 1
                                AND i.ins_start_date < '2020-01-01' AND i.ins_end_date > '2023-01-01'
                                "
) %>%
  group_by(date, birthyear) %>%
  summarise(vaccinated_total_in_day = sum(vaccinated_total_in_day, na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(age_class = assign_age_class(birthyear)) %>%
  filter(age_class != "<5")

pop_by_age_class <- dbGetQuery( con,
                                     "
                                     SELECT DISTINCT
                                     birthyear,
                                     COUNT(DISTINCT c.client_id) OVER (PARTITION BY c.birthyear) AS clients_count,
                                     'OZP' AS ins_company
                                     FROM ozp_clients c
                                     LEFT JOIN ozp_ins_start_end i ON c.client_id = i.client_id
                                     WHERE i.ins_start_date < '2020-01-01' AND i.ins_end_date > '2023-01-01'
                                     
                                     UNION
                                     
                                     SELECT DISTINCT
                                     birthyear,
                                     COUNT(DISTINCT c.client_id) OVER (PARTITION BY c.birthyear) AS clients_count,
                                     'CPZP' AS ins_company
                                     FROM cpzp_clients c
                                     LEFT JOIN cpzp_ins_start_end i ON c.client_id = i.client_id
                                     WHERE i.ins_start_date < '2020-01-01' AND i.ins_end_date > '2023-01-01'
                                     "
) %>%
  group_by(birthyear) %>%
  summarise(clients_count = sum(clients_count, na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(age_class = assign_age_class(birthyear)) %>%
  group_by(age_class) %>%
  summarise(population = sum(clients_count, na.rm = TRUE),
            .groups = "drop") %>%
  filter(age_class != "<5")

vaccination_start_date <- pop_by_age_class %>%
  select(age_class) %>%
  arrange(desc(age_class)) %>%
  mutate(vaccination_start_date =
           as.Date(
             c(
               "2021-01-15",
               "2021-03-01",
               "2021-04-14",
               "2021-04-23",
               "2021-04-28",
               "2021-05-05",
               "2021-05-10",
               "2021-05-17",
               "2021-05-24",
               "2021-05-26",
               "2021-06-04",
               "2021-07-01",
               "2021-12-13"
             )
           ))

Věkové kategorie jsou pro účely sledování proočkovanosti rozřazeny podle harmonogramu zahájení vakcinace. V některých případech byla část populace vakcinovaná dříve, pokud bylo zpřístupnění zahájeno na základě jiného kritéria než věkového (např. personál zdravotnických zařízení). Tyto věkové kategorie jsou následně použité i později při detailních analýzách.

vaccination_start_date_cz <- vaccination_start_date %>% filter(age_class != "<5")
colnames(vaccination_start_date_cz) <- c("Věková kategorie", "Datum zpřístupnění vakcinace")

vaccination_start_date_cz %>%
  kable("html", caption = "Harmonogram zahájení vakcinace") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = TRUE)
Harmonogram zahájení vakcinace
Věková kategorie Datum zpřístupnění vakcinace
80+ 2021-01-15
70-79 2021-03-01
65-69 2021-04-14
60-64 2021-04-23
55-59 2021-04-28
50-54 2021-05-05
45-49 2021-05-10
40-44 2021-05-17
35-39 2021-05-24
30-34 2021-05-26
16-29 2021-06-04
12-15 2021-07-01
5-11 2021-12-13

1.2.1 Proočkovanost populace

Je zřejmé, že každá věková kategorie se chovala v přístupu k očkování jinak. Patrný trend je v logice “Čím starší populace, tím vyšší proočkovanost”.

vaccinated_first_shot_cumulative_ratio <- vaccination_first_shot_by_birthyear %>%
  group_by(age_class, date) %>%
  summarise(
    vaccinated_in_day = sum(vaccinated_total_in_day, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(age_class, date) %>%
  group_by(age_class) %>%
  mutate(cumulative_vaccinated = cumsum(vaccinated_in_day)) %>%
  ungroup() %>%
  inner_join(pop_by_age_class, by = "age_class") %>%
  mutate(vaccinated_ratio = cumulative_vaccinated / population)

ggplot(
  vaccinated_first_shot_cumulative_ratio %>%
    filter(date <= as.Date("2022-01-01")),
  aes(x = date, y = vaccinated_ratio, color = age_class)
) +
  geom_line(size = 1) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(
    title = "Kumulativní proočkovanost dle první dávky vakcíny",
    x = "Datum",
    y = "Pročkovanost",
    color = "Věková kategorie"
  ) +
  theme_minimal()

vaccination_cumulative_relative <- vaccinated_first_shot_cumulative_ratio %>%
  left_join(vaccination_start_date, by = "age_class") %>%
  mutate(days_prior_to_vaccination_start_date = as.integer(date - vaccination_start_date))

vaccination_cumulative_relative %>%
  filter(date <= as.Date("2022-01-01")) %>%
  ggplot(aes(x = days_prior_to_vaccination_start_date, y = vaccinated_ratio, color = age_class)) +
  geom_line(size = 1) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(
    title = "Kumulativní proočkovanost",
    subtitle = "Podle doby od zahájení vakcinace pro danou kategorii",
    x = "Dnů od zahájení vakcinace",
    y = "Proočkovanost",
    color = "Věková kategorie"
  ) +
  theme_minimal(base_size = 13)

vaccinated_weekly_smooth <- vaccinated_first_shot_cumulative_ratio %>%
  filter(age_class != "5-11") %>%
  arrange(age_class, date) %>%
  group_by(age_class) %>%
  mutate(
    vaccinated_7d_avg = zoo::rollmean(
      vaccinated_in_day / population,
      k = 7,
      fill = NA,
      align = "right"
    )
  ) %>%
  ungroup() %>%
  left_join(vaccination_start_date, by = "age_class") 

ggplot(
  vaccinated_weekly_smooth %>%
    filter(date <= as.Date("2022-01-01")),
  aes(x = date, y = vaccinated_7d_avg)
) +
  geom_line(color = "steelblue", size = 0.6) +
  geom_vline(aes(xintercept = vaccination_start_date), color = "red", linetype = "dashed", size = 0.5) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  facet_wrap(~ age_class, scales = "free_y") +
  labs(
    title = "Sedmidenní průměr denní vakcinace",
    subtitle = "S počátkem vakcinace (červená)",
    x = "Datum",
    y = "% Vakcinováno"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    strip.text = element_text(face = "bold"),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

1.2.2 Eukleidovská vzdálenost proočkovanosti podle věku

Na základě míry kumulativní proočkovanosti podle doby od zahájení vakcinace lze spočítáním eukleidovských vzdáleností efektivně rozřadit věkové kategorie do širších kohort. Tyto kohorty lze později využít u analýz, které vyžadují větší množství dat (např. SCCS).

Výrazný rozdíl u skupiny 5-11 je způsoben tím, že očkování bylo zpřístupněno později a sledované období končí v roce 2022.

wide_corr_data <- vaccination_cumulative_relative %>%
  select(age_class,
         days_prior_to_vaccination_start_date,
         vaccinated_ratio) %>%
  pivot_wider(names_from = age_class, values_from = vaccinated_ratio)

mat <- wide_corr_data %>% select(-days_prior_to_vaccination_start_date) %>% as.matrix()
dist_mat <- dist(t(mat), method = "euclidean")
dist_mat_full <- as.matrix(dist_mat)

color_palette <- viridis(100)

pheatmap(
  dist_mat_full,
  color = color_palette,
  display_numbers = TRUE,
  number_color = "black",
  cluster_rows = TRUE,
  cluster_cols = TRUE,
  fontsize_number = 8,
  fontsize = 10,
  angle_col = 45,
  legend = TRUE,
  legend_breaks = seq(0, max(dist_mat_full), length.out = 5),
  legend_labels = round(seq(0, max(dist_mat_full), length.out = 5), 2),
  border_color = "grey60"
)

2 Porovnání změny trendu po roce 2020

Je nutné vzít v úvahu fakt, že samotný čas je jako prediktor obecně nedostatečný a výsledky tohoto pozorování slouží spíše k účelům nasměrování k dalšímu zkoumání.

2.1 Počet předpisů bez kategorizace

Pro předpověď vývoje dat po roce 2020 na základě předchozích dat byla použita metoda ets R knihovny forecast.

K jednoduše vysvětlitelnému propadu dochází na začátku roku 2020 a poté k nárůstu oproti očekávaným hodnotám, který nelze vysvětlit jen vyvážením předchozího lokálního extrému. Zajímavostí je, že sezónní minima jsou předpovězena s vysokou přesností a k největším odchylkám dochází u sezónních maxim.

prescriptions_clean <- prescriptions %>%
  filter(!is.na(presc_date)) %>%
  mutate(presc_date = as.Date(presc_date)) %>%
  arrange(presc_date) %>%
  filter(presc_date >= as.Date("2015-01-01") & presc_date <= as.Date("2023-12-31")) %>%
  mutate(count = 1)  


monthly_data <- prescriptions_clean %>%
  mutate(month = floor_date(presc_date, "month")) %>%
  group_by(month) %>%
  summarise(prescriptions = sum(count), .groups = 'drop')
train_data <- monthly_data %>% filter(month < as.Date("2020-01-01"))
test_data  <- monthly_data %>% filter(month >= as.Date("2020-01-01"))


ts_train <- ts(train_data$prescriptions, start = c(2015, 1), frequency = 12)


fit <- ets(ts_train)


n_ahead <- nrow(test_data)
forecast_result <- forecast(fit, h = n_ahead)


comparison <- test_data %>%
  mutate(
    predicted = forecast_result$mean,
    residual = prescriptions - predicted
  )
ggplot() +
  geom_line(data = monthly_data, aes(x = month, y = prescriptions), color = "blue", size = 1, alpha = 0.6) +
  geom_line(data = comparison, aes(x = month, y = predicted), color = "red", size = 1, linetype = "dashed") +
  geom_vline(xintercept = as.Date("2020-01-01"), linetype = "dotted") +
  labs(title = "Realita vs predikované hodnoty (ETS)",
       x = "Datum", y = "Počet předpisů")

ggplot(comparison, aes(x = month, y = residual)) +
  geom_line(color = "steelblue", size = 1) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(title = "Odchylky 2020+", x = "Datum", y = "Odchylka")

Pro otestování významnosti změny úhrnu byl proveden t-test mezi průměrnými odchylkami od ETS modelu pro data před a po roce 2020. Rozdíl je statisticky významný.

res_train <- residuals(fit)
res_test <- test_data$prescriptions - as.numeric(forecast_result$mean)

residuals_df <- bind_rows(
  tibble(month = train_data$month, residual = res_train, period = "train"),
  tibble(month = test_data$month, residual = res_test, period = "test")
)


ttest <- t.test(
  residuals_df$residual[residuals_df$period == "train"],
  residuals_df$residual[residuals_df$period == "test"]
)

cat("```\n", paste(capture.output(ttest), collapse = "\n"), "\n```")
 
    Welch Two Sample t-test

data:  residuals_df$residual[residuals_df$period == "train"] and residuals_df$residual[residuals_df$period == "test"]
t = -3.4167, df = 47, p-value = 0.001317
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -2118.0197  -548.1661
sample estimates:
    mean of x     mean of y 
-1.347016e-03  1.333092e+03 
 

2.2 Počet předpisů podle pohlaví a vakcinace

Pro analýzu užších skupin byla použita drobně odlišná metodika. V první iteraci je postup téměř identický jako v původním případu. Pro analýzu samotného trendu jsou ale data nejprve očištěna o sezónnost a šum pomocí metody stl a poté linearizována prostou lineární regresí. Tu lze zaprvé porovnat s očištěnými daty pro zjištění kvality modelu a zadruhé je změna trendu jednoduše patrná porovnáním sklonu každé z přímek.

Pro evaluaci odchylek je u neočištěných dat použito relativní RMSE (rRMSE), které je spočítáno jako RMSE po roce 2020 vydělené aritmetickým průměrem průběhu do roku 2020. To poté znázorňuje relativní odchylku od původních dat.

U druhého přístupu umožňuje linerání regrese i vyhodnocení kvality modelu. Toho je docíleno výpočtem koeficientu determinace pro data do roku 2022 (r_squared) a p-hodnoty samotného modelu (p_value_lm).

analyze_group <- function(df, group_id, group_type) {
  monthly_data <- df %>%
    mutate(month = floor_date(presc_date, "month")) %>%
    group_by(month) %>%
    summarise(prescriptions = sum(count), .groups = 'drop') %>%
    arrange(month)

  train_data <- monthly_data %>% filter(month < as.Date("2020-01-01"))
  test_data  <- monthly_data %>% filter(month >= as.Date("2020-01-01"))

  if (nrow(train_data) < 24 || nrow(test_data) < 12) {
    return(NULL)
  }

  train_records <- df %>% filter(presc_date < as.Date("2020-01-01")) %>% nrow()

  ts_train <- ts(train_data$prescriptions, start = c(year(min(train_data$month)), month(min(train_data$month))), frequency = 12)

 fit <- tryCatch(ets(ts_train, model = "ZZM"), error = function(e) NULL)
  if (is.null(fit)) return(NULL)

  forecast_result <- forecast(fit, h = nrow(test_data))
  residuals_train <- residuals(fit)
  residuals_test <- test_data$prescriptions - as.numeric(forecast_result$mean)

  ttest <- t.test(residuals_train, residuals_test)
  rmse_train <- sqrt(mean(residuals_train^2))
  rmse <- sqrt(mean(residuals_test^2))
  p_train <- residuals(fit)
  
  mean_train <- mean(train_data$prescriptions)
  rrmse <- rmse / mean_train

  plot_data <- bind_rows(
    train_data %>% mutate(predicted = as.numeric(fitted(fit))),
    test_data %>% mutate(predicted = as.numeric(forecast_result$mean))
  ) %>%
    mutate(group = group_id, group_type = group_type)

  tibble(
    group = group_id,
    group_type = group_type,
    train_records = train_records,
    rmse_train = rmse_train,
    rmse = rmse,
    rrmse = rrmse,
    t_stat = ttest$statistic,
    p_value = ttest$p.value,
    mean_train_resid = mean(residuals_train),
    mean_test_resid = mean(residuals_test),
    plot_data = list(plot_data)
  )
}

analyze_group_trend <- function(df, group_id, group_type) {
  monthly_data <- df %>%
    mutate(month = floor_date(presc_date, "month")) %>%
    group_by(month) %>%
    summarise(prescriptions = sum(count), .groups = "drop") %>%
    arrange(month)

  if (nrow(monthly_data) < 36) return(NULL)

  ts_full <- ts(monthly_data$prescriptions,
                start = c(year(min(monthly_data$month)), month(min(monthly_data$month))),
                frequency = 12)
  
  stl_fit <- stl(ts_full, s.window = "periodic", robust = TRUE)
  trend_full <- stl_fit$time.series[, "trend"]

  trend_df <- tibble(
    month = monthly_data$month,
    trend = as.numeric(trend_full)
  ) %>%
    mutate(
      time_index = 1:n(),
      post_2020 = month >= as.Date("2020-01-01")
    )
  
  pre_df <- filter(trend_df, !post_2020)
  post_df <- filter(trend_df, post_2020)


  lm_pre <- lm(trend ~ time_index, data = pre_df)
  lm_post <- lm(trend ~ time_index, data = post_df)


  lm_pre_summary <- summary(lm_pre)
  pre_p_value <- lm_pre_summary$coefficients[2, 4]
  pre_r_squared <- lm_pre_summary$r.squared


  post_df <- post_df %>%
    mutate(
      extrapolated_pre_trend = predict(lm_pre, newdata = post_df),
      actual_post_trend = predict(lm_post, newdata = post_df),
      residuals = actual_post_trend - extrapolated_pre_trend
    )


  rmse_post <- sqrt(mean(post_df$residuals^2))
  mean_pre_trend <- mean(pre_df$trend)
  rrmse_post <- rmse_post / mean_pre_trend

  ttest <- t.test(post_df$residuals)

  trend_df <- trend_df %>%
    mutate(
      pre_fit = predict(lm_pre, newdata = trend_df),
      post_fit = if_else(post_2020, predict(lm_post, newdata = trend_df), NA_real_)
    )

  tibble(
    group = group_id,
    group_type = group_type,
    n_records = nrow(df),
    slope_pre = coef(lm_pre)["time_index"],
    slope_post = coef(lm_post)["time_index"],
    p_value_pre = pre_p_value,
    r_squared_pre = pre_r_squared,
    rmse_post = rmse_post,
    rrmse_post = rrmse_post,
    mean_residual = mean(post_df$residuals),
    t_stat = ttest$statistic,
    p_value = ttest$p.value,
    plot_data = list(trend_df)
  )
}
results_vacc <- prescriptions_clean %>%
  filter(!is.na(vaccinated)) %>%
  group_by(vaccinated) %>%
  group_split() %>%
  map(~ {
    group_val <- unique(.x$vaccinated)
    analyze_group(df = filter(prescriptions_clean, vaccinated == group_val), 
                  group_id = group_val, 
                  group_type = "Vaccination")
  }) %>%
  compact() %>%
  bind_rows()

results_sex <- prescriptions_clean %>%
  filter(!is.na(sex)) %>%
  group_by(sex) %>%
  group_split() %>%
  map(~ {
    group_val <- unique(.x$sex)
    analyze_group(df = filter(prescriptions_clean, sex == group_val), 
                  group_id = group_val, 
                  group_type = "Sex")
  }) %>%
  compact() %>%
  bind_rows()

results_all <- bind_rows(results_vacc, results_sex) %>%
  arrange(p_value)

plot_data_all <- results_all %>%
  unnest(plot_data, names_sep = "_")

ETS model vykazuje největší rozdíl u nevakcinované skupiny a u mužů. T-test rozdílu odchylek je v obou případech statisticky významný.

results_all %>%
  mutate(Skupina = group, RMSE = rmse, rRMSE = rrmse) %>%
  select(
    Skupina, RMSE, rRMSE, t_stat, p_value
  ) %>%
  kable(format = "html", caption = "Porovnání předpisovosti v čase (ETS)") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání předpisovosti v čase (ETS)
Skupina RMSE rRMSE t_stat p_value
M 1396.8746 0.1692402 -5.4883544 0.0000016
Bez vakcinace 696.6278 0.1651938 -4.4321352 0.0000557
F 1759.1357 0.1127051 -1.4953911 0.1414972
Vakcinace 2068.4206 0.1052896 0.0251487 0.9800428
ggplot(plot_data_all, aes(x = plot_data_month)) +
  geom_line(aes(y = plot_data_prescriptions), color = "blue") +
  geom_line(aes(y = plot_data_predicted), color = "red", linetype = "dashed") +
  geom_vline(xintercept = as.Date("2020-01-01"), linetype = "dotted") +
  facet_wrap(~ paste(plot_data_group_type, plot_data_group, sep = ": "), scales = "free_y") +
  labs(
    title = "Realita vs Předpověď (ETS)",
    subtitle = "Modrá = Realita, Červená = Předpověď",
    x = "Datum", y = "Počet předpisů"
  ) +
  theme_minimal()

Linearizovaný trend vykazuje nárůst od předpokladu ve všech kategoriích. Největší rozdíl je opět patrný u skupiny neočkovaných a mužů. U nevakcinované populace je nutné vzít v úvahu velmi nízký koeficient determinace a tedy i nízké prediktivní schopnosti samotného modelu.

results_vacc <- prescriptions_clean %>%
  filter(!is.na(vaccinated)) %>%
  group_by(vaccinated) %>%
  group_split() %>%
  map(~ analyze_group_trend(.x, group_id = unique(.x$vaccinated), group_type = "Vaccinated")) %>%
  compact() %>%
  bind_rows()

results_sex <- prescriptions_clean %>%
  filter(!is.na(sex)) %>%
  group_by(sex) %>%
  group_split() %>%
  map(~ analyze_group_trend(.x, group_id = unique(.x$sex), group_type = "Sex")) %>%
  compact() %>%
  bind_rows()

results_all <- bind_rows(results_vacc, results_sex)

plot_data_all <- results_all %>%
  unnest(plot_data)
results_all %>%
  mutate(Skupina = group, RMSE = rmse_post, rRMSE = rrmse_post, r_squared = r_squared_pre, p_value_diff = p_value, p_value_lm = p_value_pre) %>%
  select(
    Skupina, RMSE, rRMSE, p_value_lm, r_squared
  ) %>%
  kable(format = "html", caption = "Porovnání předpisovosti v čase (linearizovaný trend)") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání předpisovosti v čase (linearizovaný trend)
Skupina RMSE rRMSE p_value_lm r_squared
Bez vakcinace 565.2184 0.1337668 0.0004409 0.1932531
Vakcinace 1042.0182 0.0528575 0.0000000 0.9306673
F 860.4047 0.0549194 0.0000000 0.9032496
M 765.4342 0.0924169 0.0000000 0.8175001
ggplot(plot_data_all, aes(x = month)) +
  geom_line(aes(y = trend), color = "black", size = 0.8) +               
  geom_line(aes(y = pre_fit), color = "blue", size = 1.2, linetype = "dashed") +   
  geom_line(aes(y = post_fit), color = "red", size = 1.2, linetype = "dashed") +   
  geom_vline(xintercept = as.Date("2020-01-01"), linetype = "solid", color = "gray40") +
  facet_wrap(~ paste(group_type, group, sep = ": "), scales = "free_y") +
  labs(
    title = "Linearizované trendy v předpisovosti",
    subtitle = "Modrá = Linearizovaný trend do 2020; Červená = Linearizovaný trend po 2020",
    x = "Month",
    y = "Počet předpisů"
  ) +
  theme_minimal() +
  theme(
    strip.text = element_text(face = "bold"),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

2.3 Počet předpisů podle věkové kategorie

Pro porovnání věkových kategorií byla použita stejná metodika jako u skupin podle vakcinace a pohlaví.

results_age <- prescriptions_clean %>%
  filter(!is.na(age_class)) %>%
  group_by(age_class) %>%
  group_split() %>%
  map(~ {
    group_val <- unique(.x$age_class)
    analyze_group(df = filter(prescriptions_clean, age_class == group_val), 
                  group_id = group_val, 
                  group_type = "Age Class")
  }) %>%
  compact() %>%
  bind_rows()

results_age <- results_age %>%
  arrange(group)

plot_data_age <- results_age %>%
  unnest(plot_data, names_sep = "_")

K největšímu relativnímu rozdílu dochází u skupin 5-11, 12-15 a 16-29. To částečně vysvětluje i předchozí porovnání, kde největší rozdíl vykazovala skupina nevakcinovaných. Rozdíl reziduál je však významný jen u skupiny 16-29. Skupiny 60-64 a 65-69 vykazují také signifikantní odchylku od očekávaných hodnot.

results_age %>%
  mutate(Skupina = group, RMSE = rmse, rRMSE = rrmse) %>%
  select(
    Skupina, RMSE, rRMSE, t_stat, p_value
  ) %>%
  kable(format = "html", caption = "Porovnání předpisovosti v čase (ETS)") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání předpisovosti v čase (ETS)
Skupina RMSE rRMSE t_stat p_value
5-11 116.56860 0.2761308 0.1210174 0.9041932
12-15 60.02150 0.2807804 -1.1662273 0.2494064
16-29 238.75694 0.2524881 -4.8869924 0.0000123
30-34 87.74207 0.1384490 7.4119505 0.0000000
35-39 106.79332 0.1268053 -2.5359254 0.0145986
40-44 208.35234 0.1473983 -4.2292198 0.0001075
45-49 290.60068 0.1308266 -4.2824696 0.0000905
50-54 242.86926 0.1050056 0.3906080 0.6978519
55-59 262.05711 0.0990215 0.6870654 0.4954179
60-64 312.48431 0.1194154 -4.6045264 0.0000316
65-69 495.21694 0.1727451 -4.7074821 0.0000224
70-79 618.10440 0.1224329 0.1086331 0.9139560
80+ 253.20135 0.1508588 5.8897066 0.0000004
ggplot(plot_data_age, aes(x = plot_data_month)) +
  geom_line(aes(y = plot_data_prescriptions), color = "blue") +
  geom_line(aes(y = plot_data_predicted), color = "red", linetype = "dashed") +
  geom_vline(xintercept = as.Date("2020-01-01"), linetype = "dotted") +
  facet_wrap(~ plot_data_group, scales = "free_y", ncol = 2) +
  labs(
    title = "Realita vs Předpověď (ETS)",
    subtitle = "Modrá = Realita, Červená = Předpověď",
    x = "Datum", y = "Počet předpisů"
  ) +
  theme_minimal()

results_age <- prescriptions_clean %>%
  filter(!is.na(age_class)) %>%
  group_by(age_class) %>%
  group_split() %>%
  map(~ {
    group_val <- unique(.x$age_class)
    analyze_group_trend(df = filter(prescriptions_clean, age_class == group_val), 
                  group_id = group_val, 
                  group_type = "Age Class")
  }) %>%
  compact() %>%
  bind_rows()

results_age <- results_age %>%
  arrange(group)

plot_data_age <- results_age %>%
  unnest(plot_data)

Porovnání linearizovaných trendů víceméně potvrzuje některé předpoklady nalezené v předchozí analýze. Bohužel mladší skupiny (pod 15 let) mají natolik výrazné výkyvy i po očištění o sezónnost, že není možné jejich průběh linearizovat s vysokým koeficientem determinace. Skupina 16-29 již nevykazuje tak výraznou odchylku od předpokladu a u starší populace vykazuje výrazný rozdíl jen 65-69.

results_age %>%
  mutate(Skupina = group, RMSE = rmse_post, rRMSE = rrmse_post, r_squared = r_squared_pre, p_value_diff = p_value, p_value_lm = p_value_pre) %>%
  select(
    Skupina, RMSE, rRMSE, p_value_lm, r_squared
  ) %>%
  kable(format = "html", caption = "Porovnání předpisovosti v čase (linearizovaný trend)") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání předpisovosti v čase (linearizovaný trend)
Skupina RMSE rRMSE p_value_lm r_squared
5-11 223.26393 0.5500878 0.0000000 0.6525649
12-15 95.46844 0.4596474 0.0000000 0.6818070
16-29 65.26420 0.0689561 0.0000000 0.9517284
30-34 41.46398 0.0652822 0.0000000 0.9574665
35-39 43.78888 0.0518640 0.0000000 0.9903928
40-44 105.50875 0.0746866 0.0000000 0.9401286
45-49 239.89668 0.1077107 0.0000000 0.9644207
50-54 155.70665 0.0672456 0.0000000 0.9816383
55-59 130.29275 0.0491745 0.0000000 0.9537575
60-64 117.37329 0.0447723 0.0000000 0.8245083
65-69 341.38374 0.1183617 0.0000000 0.5065845
70-79 192.17282 0.0379005 0.0000000 0.8453986
80+ 103.83404 0.0616256 0.1008976 0.0457235
ggplot(plot_data_age, aes(x = month)) +
  geom_line(aes(y = trend), color = "black", size = 0.8) +               
  geom_line(aes(y = pre_fit), color = "blue", size = 1.2, linetype = "dashed") +   
  geom_line(aes(y = post_fit), color = "red", size = 1.2, linetype = "dashed") +   
  geom_vline(xintercept = as.Date("2020-01-01"), linetype = "solid", color = "gray40") +
  facet_wrap(~ group, scales = "free_y", ncol = 2) +
  labs(
    title = "Linearizované trendy v předpisovosti",
    subtitle = "Modrá = Linearizovaný trend do 2020; Červená = Linearizovaný trend po 2020",
    x = "Month",
    y = "Počet předpisů"
  ) +
  theme_minimal() +
  theme(
    strip.text = element_text(face = "bold"),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

Tato úvodní analýza poskytla vhodný startovací bod pro následné zkoumání konkrétních věkových kohort. Konkrétní skupiny vybrané pro další zkoumání tedy jsou: 12-15 z důvodu největšího výkyvu (skupina je ale bohužel natolik malá, že pravděpodobně nebude možné dojít k relevantním závěrům), 16-29, 30-39 a 40-59. Skupina 5-11 je vyřazena z několika důvodů. Prvním je nižší proočkovanost a druhým je fakt, že nejmladší populace této skupiny je na začátku dlouhodobého meření buď nenarozená nebo čerstvě na světě a může zde probíhat mnoho výkyvů. Skupina 60-69 by jistě také poskytla zajímavé výsledky, ale z časových důvodů a z důvodu zadání soutěže (orientace především na mladší populaci) byla prozatím také vyřazena.

3 Metodika detailní analýzy

Detailní analýza je rozdělená na dvě části. Dlouhodobou a krátkodobou. Zatímco dlouhodobá analýza porovnává trendy v období před covidem (<2020) vs po covidu (>=2022), tak krátkodobá analýza sleduje efekty v horizontu jednoho roku od vakcíny s porovnáním stejně dlouhého období před vakcinací.

Pokud dochází k výkyvu celkové předpisovosti, tak musí zákonitě docházet i ve výkyvu buď v prvopředpisovosti, multipředpisovosti nebo obojím. Tyto dvě metriky jsou v dlouhodobé analýze měřené na 1000 osob v riziku, aby byly alespoň částečně porovnatelné. Pro prvopředpisovost tedy počet prvopředpisů na 1000 dosud nepředepsaných osob a pro multipředpisovost poměr multipředpisovaných osob v daném období na 1000 předpisovaných osob (období jsou roční a měsíční). Zároveň je také sledovaná míra retence (tedy kolik procent prvopředepisovaných bylo předepsaných i v následujícíh letech). Z těchto naměřených hodnot je následně vypočten poměr po vs před (a t-test na rozdíl průměrných hodnot) a statistická významnost rozdílu mezi očkovanými a neočkovanými je vyhodnocena pomocí p-hodnoty pro postCOVID:vaccinated lineárního modelu log(RatePer1000) ~ postCOVID * vaccinated, kde proměnné postCOVID a vaccinated nabývají hodnot pouze 1 a 0.

Pro krátkodobé efekty je použita metoda SCCS, konkrétně tedy standardsccs(). Očkovaní a neočkovaní jsou zde porovnáni s vlastní baseline 365 dní před vakcinací a následně jsou posouzena rizika prvopředpisů a předpisů celkem po vakcíně. Rizika jsou následně porovnána mezi očkovanými a neočkovanými. V případě neočkovaných je datum vakcinace nastaveno jako rozhodné datum v dané skupině. Skupina očkovaných vybraná pro analýzu je podle vakcinace (první dávka) +- 7 dní od rozhodného data. Do analýzy jsou oproti všem ostatním přístupům zahrnuti klienti, kteří byli pojištěni po dobu +- jednoho roku od vakcinace (v ostatních přístupech jsou zahrnuti klienti pojištěni minimálně po celých 8 let) a multipředpisy injekcí / infuzí v rozsahu 10 dnů jsou evidovány jako pouze jeden předpis, aby měla každá událost zhruba stejnou váhu. Období po podání vakcíny byla nastavena jako 0-29, 30-89, 90-179, 180-269 a 270-365 dní. Kratší intervaly bohužel často vykazovaly nižší množství dat a tedy insignifikantní výsledky. Skupiny 16-29 a 30-39 zde byly sloučeny do jedné z důvodu nízké relevance výstupů u samostatné evaluace. Pro testování rozdílů mezi skupinami byl použit Waldův test na rozdíl log(Hazard Ratio) v rámci DiD analýzy.

age_cohorts <- vaccination_start_date %>%
  arrange(age_class)

age_cohorts <- age_cohorts[2:9 , ]

age_cohorts <- age_cohorts %>%
  mutate(
      cohort = c(rep("12-15", 1), rep("16-39", 3), rep("40-59", 4))
  )


vaccination_peaks <- vaccinated_first_shot_cumulative_ratio %>%
  inner_join(age_cohorts, by = "age_class") %>%
  group_by(cohort, date) %>%
  summarise(vaccinated_in_day = sum(vaccinated_in_day, na.rm = TRUE), .groups = "drop") %>%
  arrange(cohort, date) %>%
  group_by(cohort) %>%
  mutate(vaccinated_in_day_ma7 = zoo::rollmean(vaccinated_in_day, k = 7, fill = NA, align = "right")) %>%
  ungroup()

peak_dates <- vaccination_peaks %>%
  group_by(cohort) %>%
  filter(vaccinated_in_day == max(vaccinated_in_day, na.rm = TRUE)) %>%
  slice(1) %>% 
  ungroup() %>%
  select(cohort, peak_date = date) %>%
  mutate(decisive_date = peak_date + days(7))

Z důvodu vrcholu očkování krátce po zahájení je rozhodné datum nastavené o 7 dní později než datum vrcholu.

age_cohorts %>%
  group_by(cohort) %>%
  summarise(vaccination_start_date = max(vaccination_start_date)) %>%
  left_join(peak_dates, by = "cohort") %>%
  transmute(Kohorta = cohort, `Datum zahájení vakcinace` = vaccination_start_date, `Vrchol vakcinace` = peak_date, `Rozhodné datum` = decisive_date) %>%
  kable(format = "html", caption = "Zahájení a vrcholy vakcinace dle kohorty") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Zahájení a vrcholy vakcinace dle kohorty
Kohorta Datum zahájení vakcinace Vrchol vakcinace Rozhodné datum
12-15 2021-07-01 2021-07-28 2021-08-04
16-39 2021-06-04 2021-06-08 2021-06-15
40-59 2021-05-17 2021-05-14 2021-05-21
ggplot(vaccination_peaks, aes(x = date, y = vaccinated_in_day_ma7)) +
  geom_line(color = "steelblue") +
  facet_wrap(~ cohort, scales = "free_y") +
  geom_vline(data = peak_dates, aes(xintercept = as.numeric(decisive_date)), color = "red", linetype = "dashed") +
  labs(title = "Sedmidenní průměr denní vakcinace a vrcholy vakcinace",
       x = "Datum", y = "Vakcinovaných za den") +
  theme_minimal()

4 Věková kategorie 12-15

4.1 Dlouhodobý efekt

4.1.1 Distribuce věku a pohlaví

Ačkoliv je věková distribuce podle vakcinace mírně odlišná, neměla by mít výrazný vliv z pohledu rozdílných rizik podle stáří.

unique_clients <- prescriptions %>%
  filter(!is.na(age_2021), !is.na(sex), !is.na(vaccinated)) %>%
  distinct(client_id, .keep_all = TRUE)  

age_data <- unique_clients %>%
  group_by(vaccinated, age_2021) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(vaccinated) %>%
  mutate(prop = count / sum(count),
         variable = "Age",
         value = as.character(age_2021)) %>%
  ungroup()

sex_data <- unique_clients %>%
  group_by(vaccinated, sex) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(vaccinated) %>%
  mutate(prop = count / sum(count),
         variable = "Sex",
         value = sex) %>%
  ungroup()
ggplot(unique_clients, aes(x = age_2021, fill = factor(vaccinated))) +
  geom_histogram(aes(y = after_stat(density)), position = "identity", alpha = 0.4, binwidth = 1) +
  scale_fill_brewer(palette = "Set1", name = "Vaccinated") +
  scale_x_continuous(breaks = seq(floor(min(unique_clients$age)), ceiling(max(unique_clients$age)), by = 1)) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    title = "Věkové rozložení podle statusu vakcinace",
    x = "Věk", y = "%",
    fill = "Vakcinace"
  ) +
  theme_minimal()

ggplot(sex_data, aes(x = sex, y = prop, color = vaccinated, group = vaccinated)) +
  geom_line(size = 1) +
  geom_point(size = 2) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    title = "Pohlaví podle vakcinace",
    x = "Pohlaví",
    y = "%",
    color = "Vakcinace"
  ) +
  theme_minimal()

4.1.2 Předpověď vs realita

ETS model vyhodnocuje rozdíl v nastaveném trendu zejména u ženské a neočkované populace. Tyto rozdíly po vs před rokem 2020 jsou statisticky významné.

prescriptions_clean <- prescriptions %>%
  filter(!is.na(presc_date)) %>%
  mutate(presc_date = as.Date(presc_date)) %>%
  arrange(presc_date) %>%
  filter(presc_date >= as.Date("2015-01-01") & presc_date <= as.Date("2023-12-31")) %>%
  mutate(count = 1)  

results_vacc <- prescriptions_clean %>%
  filter(!is.na(vaccinated)) %>%
  group_by(vaccinated) %>%
  group_split() %>%
  map(~ {
    group_val <- unique(.x$vaccinated)
    analyze_group(df = filter(prescriptions_clean, vaccinated == group_val), 
                  group_id = group_val, 
                  group_type = "Vaccination")
  }) %>%
  compact() %>%
  bind_rows()

results_sex <- prescriptions_clean %>%
  filter(!is.na(sex)) %>%
  group_by(sex) %>%
  group_split() %>%
  map(~ {
    group_val <- unique(.x$sex)
    analyze_group(df = filter(prescriptions_clean, sex == group_val), 
                  group_id = group_val, 
                  group_type = "Sex")
  }) %>%
  compact() %>%
  bind_rows()

results_all <- bind_rows(results_vacc, results_sex) %>%
  arrange(p_value)

plot_data_all <- results_all %>%
  unnest(plot_data, names_sep = "_")
results_all %>%
  mutate(Skupina = group, RMSE = rmse, rRMSE = rrmse) %>%
  select(
    Skupina, RMSE, rRMSE, t_stat, p_value
  ) %>%
  kable(format = "html", caption = "Porovnání předpisovosti v čase (ETS)") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání předpisovosti v čase (ETS)
Skupina RMSE rRMSE t_stat p_value
Bez vakcinace 47.67759 0.4239264 -5.4861625 0.0000016
F 36.48232 0.4145718 -3.1989369 0.0024708
M 32.17169 0.2558046 1.5519815 0.1273749
Vakcinace 33.29083 0.3286361 -0.8040151 0.4254362
ggplot(plot_data_all, aes(x = plot_data_month)) +
  geom_line(aes(y = plot_data_prescriptions), color = "blue") +
  geom_line(aes(y = plot_data_predicted), color = "red", linetype = "dashed") +
  geom_vline(xintercept = as.Date("2020-01-01"), linetype = "dotted") +
  facet_wrap(~ paste(plot_data_group_type, plot_data_group, sep = ": "), scales = "free_y") +
  labs(
    title = "Realita vs Předpověď (ETS)",
    subtitle = "Modrá = Realita, Červená = Předpověď",
    x = "Datum", y = "Počet předpisů"
  ) +
  theme_minimal()

4.1.3 Prvopředpisovost

Všeobecná prvopředpisovost je na celé populaci mírně vyšší v období po covidu.

first_presc <- dbGetQuery(con,
                          "
                          SELECT
                            c.client_id,
                            c.n_vacc,
                            c.n_presc,
                            c.sex,
                            c.birthyear,
                            FIRST(p.date) OVER (PARTITION BY c.client_id ORDER BY p.date ASC) AS first_presc_date,
                            'CPZP' AS ins_company
                          FROM
                            cpzp_clients c
                          LEFT JOIN cpzp_ins_start_end i ON c.client_id = i.client_id
                          LEFT JOIN cpzp_prescriptions p ON c.client_id = p.client_id AND p.presc_order = 1
                          WHERE
                            i.ins_start_date < '2015-01-01' AND i.ins_end_date = '9999-12-31'
                            AND 2021 - c.birthyear BETWEEN 12 AND 15
                            
                          UNION
                          
                          SELECT
                            c.client_id,
                            c.n_vacc,
                            c.n_presc,
                            c.sex,
                            c.birthyear,
                            FIRST(p.date) OVER (PARTITION BY c.client_id ORDER BY p.date ASC) AS first_presc_date,
                            'OZP' AS ins_company
                          FROM
                            ozp_clients c
                          LEFT JOIN ozp_ins_start_end i ON c.client_id = i.client_id
                          LEFT JOIN ozp_prescriptions p ON c.client_id = p.client_id AND p.presc_order = 1
                          WHERE
                            i.ins_start_date < '2015-01-01' AND i.ins_end_date IS NULL
                            AND 2021 - c.birthyear BETWEEN 12 AND 15
                          ") %>%
  mutate(first_presc_yearmonth = floor_date(first_presc_date, unit = "month"),
         sex = ifelse(sex == "Z", "F", sex),
         vacc_status = ifelse(n_vacc > 0, "Vakcinace", "Bez vakcinace"),
         age_2021 = 2021 - birthyear)

monthly_first_presc <- first_presc %>%
  count(first_presc_yearmonth, name = "new_prescriptions") %>%
  arrange(first_presc_yearmonth) %>%
  mutate(
    cum_new_presc = cumsum(new_prescriptions),  
    total_population = nrow(first_presc),
    at_risk_population = total_population - lag(cum_new_presc, default = 0),
    rate = new_prescriptions / at_risk_population
  ) %>%
  mutate(
    period_group = case_when(
      first_presc_yearmonth < as.Date("2020-01-01") ~ "Pre-COVID",
      first_presc_yearmonth > as.Date("2021-12-31") ~ "Post-COVID",
      TRUE ~ "COVID-Shock"
    )
  )

monthly_clean <- monthly_first_presc %>%
  filter(first_presc_yearmonth >= as.Date("2018-01-01")) %>%
  filter(period_group %in% c("Pre-COVID", "Post-COVID"))

ttest <- t.test(rate ~ period_group, data = monthly_clean)
ggplot(monthly_first_presc %>% filter(first_presc_yearmonth >= as.Date("2018-01-01")),
       aes(x = first_presc_yearmonth, y = rate * 1000)) +
  geom_line(color = "steelblue", size = 1) +
  annotate("rect",
           xmin = as.Date("2020-01-01"), xmax = as.Date("2021-12-31"),
           ymin = -Inf, ymax = Inf, alpha = 0.2, fill = "gray50") +
  scale_y_continuous(name = "Počet prvopředpisů na 1000 nepředepsaných klientů", limits = c(0, 2)) +
  labs(
    title = "Prvopředpisovost v čase",
    subtitle = "Šedá = covidové období (vyřazeno z t-testu)",
    x = "Datum"
  ) +
  theme_minimal()

cat("```\n", paste(capture.output(ttest), collapse = "\n"), "\n```")
 
    Welch Two Sample t-test

data:  rate by period_group
t = 2.9636, df = 44.846, p-value = 0.004853
alternative hypothesis: true difference in means between group Post-COVID and group Pre-COVID is not equal to 0
95 percent confidence interval:
 5.783436e-05 3.032637e-04
sample estimates:
mean in group Post-COVID  mean in group Pre-COVID 
             0.001248165              0.001067616 
 
total_by_vacc <- first_presc %>%
  count(vacc_status, name = "total_population")

monthly_by_vacc <- first_presc %>%
  count(vacc_status, first_presc_yearmonth, name = "new_prescriptions") %>%
  arrange(vacc_status, first_presc_yearmonth) %>%
  group_by(vacc_status) %>%
  mutate(
    cum_prescriptions = cumsum(new_prescriptions),
    total_population = total_by_vacc$total_population[match(vacc_status, total_by_vacc$vacc_status)],
    at_risk_population = total_population - lag(cum_prescriptions, default = 0),
    rate = new_prescriptions / at_risk_population
  ) %>%
  ungroup() %>%
  mutate(
    period_group = case_when(
      first_presc_yearmonth < as.Date("2020-01-01") ~ "Pre-COVID",
      first_presc_yearmonth > as.Date("2021-12-31") ~ "Post-COVID",
      TRUE ~ "COVID-Shock"
    )
  )
  
monthly_by_vacc_clean <- monthly_by_vacc %>%
  filter(first_presc_yearmonth >= as.Date("2018-01-01")) %>%
  filter(period_group %in% c("Pre-COVID", "Post-COVID"))

first_presc_table <- monthly_by_vacc_clean %>%
  group_by(vacc_status) %>%
  summarise(
    p_value = t.test(rate ~ period_group)$p.value,
    mean_before = mean(rate[period_group == "Pre-COVID"]),
    mean_after = mean(rate[period_group == "Post-COVID"]),
    rate = mean_after / mean_before
  )

Prvopředpisovost po covidu vykazuje nárůst u obou kategorií, v případě očkovaných však vyšší a statisticky významný.

first_presc_table %>%
  transmute(Vakcinace = vacc_status, p_value, `Průměr před` = mean_before * 1000, `Průměr po` = mean_after * 1000, `Poměr` = rate) %>%
  kable(format = "html", caption = "Porovnání prvopředpisovosti v čase podle vakcinace") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání prvopředpisovosti v čase podle vakcinace
Vakcinace p_value Průměr před Průměr po Poměr
Bez vakcinace 0.1763358 1.010252 1.094448 1.083342
Vakcinace 0.0007757 1.132464 1.423175 1.256706
ggplot(monthly_by_vacc %>% filter(first_presc_yearmonth >= as.Date("2018-01-01")),
       aes(x = first_presc_yearmonth, y = rate * 1000, color = vacc_status)) +
  geom_line(size = 1) +
  annotate("rect",
           xmin = as.Date("2020-01-01"), xmax = as.Date("2021-12-31"),
           ymin = -Inf, ymax = Inf, alpha = 0.2, fill = "gray50") +
  scale_y_continuous(name = "Prescriptions per 1,000 at-risk clients", limits = c(0, 2.5)) +
  labs(
    title = "Prvopředpisovost v čase podle vakcinace",
    subtitle = "Šedá = covidové období (vyřazeno z t-testu)",
    x = "Datum",
    color = "Vakcinace"
  ) +
  theme_minimal()

DiD model ovšem statisticky významný rozdíl mezi skupinami nepotvrzuje (p hodnota pro interakci vaccinated:post_covid).

monthly_rates <- monthly_by_vacc %>%
  filter(period_group %in% c("Pre-COVID", "Post-COVID")) %>%
  mutate(
    vaccinated = if_else(vacc_status == "Vakcinace", 1, 0),
    post_covid = if_else(first_presc_yearmonth >= as.Date("2022-01-01"), 1, 0),
    interaction = vaccinated * post_covid)

did_model <- lm(log(rate * 1000) ~ vaccinated * post_covid, data = monthly_rates)

cat("```\n", paste(capture.output(summary(did_model)), collapse = "\n"), "\n```")
 
Call:
lm(formula = log(rate * 1000) ~ vaccinated * post_covid, data = monthly_rates)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.80598 -0.26332  0.00375  0.17636  1.17803 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)            0.22949    0.04683   4.900 2.28e-06 ***
vaccinated             0.09257    0.06623   1.398   0.1641    
post_covid            -0.15655    0.08762  -1.787   0.0758 .  
vaccinated:post_covid  0.16884    0.12391   1.363   0.1749    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3628 on 164 degrees of freedom
Multiple R-squared:  0.05497,   Adjusted R-squared:  0.03768 
F-statistic:  3.18 on 3 and 164 DF,  p-value: 0.02555
 

4.1.4 Multipředpisovost

4.1.4.1 Míra retence (meziroční multipředpisovost)

presc_id_counts <- dbGetQuery(con,
                          "
                          SELECT
                            c.client_id,
                            c.n_vacc,
                            c.n_presc,
                            c.sex,
                            c.birthyear,
                            2021 - c.birthyear AS age_2021,
                            p.date as presc_date,
                            p.presc_order,
                            'CPZP' AS ins_company,
                            YEAR(FIRST(p.date) OVER (PARTITION BY c.client_id ORDER BY p.presc_order ASC)) AS first_presc_year
                          FROM
                            cpzp_clients c
                          LEFT JOIN cpzp_ins_start_end i ON c.client_id = i.client_id
                          LEFT JOIN cpzp_prescriptions p ON c.client_id = p.client_id
                          WHERE
                            i.ins_start_date < '2015-01-01' AND i.ins_end_date = '9999-12-31'
                            AND 2021 - c.birthyear BETWEEN 12 AND 15
                            
                          UNION
                          
                          SELECT
                            c.client_id,
                            c.n_vacc,
                            c.n_presc,
                            c.sex,
                            c.birthyear,
                            2021 - c.birthyear AS age_2021,
                            p.date as presc_date,
                            p.presc_order,
                            'OZP' AS ins_company,
                            YEAR(FIRST(p.date) OVER (PARTITION BY c.client_id ORDER BY p.presc_order ASC)) AS first_presc_year
                          FROM
                            ozp_clients c
                          LEFT JOIN ozp_ins_start_end i ON c.client_id = i.client_id
                          LEFT JOIN ozp_prescriptions p ON c.client_id = p.client_id
                          WHERE
                            i.ins_start_date < '2015-01-01' AND i.ins_end_date IS NULL
                            AND 2021 - c.birthyear BETWEEN 12 AND 15
                          ")%>%
  mutate(presc_year = as.integer(format(as.Date(presc_date), "%Y")),
         vacc_status = ifelse(n_vacc > 0, "Vakcinace", "Bez vakcinace"))

monthly_unique_clients <- presc_id_counts %>%
  mutate(yearmonth = floor_date(presc_date, unit = "month")) %>%
  distinct(client_id, yearmonth, vacc_status) %>%
  count(yearmonth, vacc_status, name = "n_unique_clients")

presc_years <- presc_id_counts %>%
  mutate(
    presc_year = year(presc_date)
  ) %>%
  distinct(client_id, first_presc_year, presc_year, vacc_status) %>%
  filter(first_presc_year >= 2017, first_presc_year <= 2023)


retention_data <- presc_years %>%
  mutate(year_offset = presc_year - first_presc_year) %>%
  filter(year_offset >= 0)


retention_summary <- retention_data %>%
  group_by(first_presc_year, year_offset, vacc_status) %>%
  summarise(n_clients = n_distinct(client_id), .groups = "drop")


cohort_sizes <- retention_data %>%
  filter(year_offset == 0) %>%
  group_by(first_presc_year, vacc_status) %>%
  summarise(cohort_size = n_distinct(client_id), .groups = "drop")


retention_summary <- retention_summary %>%
  left_join(cohort_sizes, by = c("first_presc_year", "vacc_status")) %>%
  mutate(retention_rate = n_clients / cohort_size)

retention_summary_filtered <- retention_summary %>%
  filter(year_offset > 0)

Vzhledem k tomu, že roční retence lze sledovat jen v několikaletém horizontu, nemáme dostatek relevantních dat pro posouzení vývoje po vs před covidem. Všeobecně je ale míra retence srovnatelná.

ggplot(retention_summary_filtered, aes(x = year_offset, y = retention_rate, color = as.factor(first_presc_year))) +
  geom_line(size = 1.2) +
  geom_point(size = 2) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_color_viridis_d(name = "Rok prvopředpisu") +
  labs(
    title = "Retence (meziroční předpisovost) podle vakcinace",
    x = "Počet let od prvopředpisu",
    y = "Míra multipředpisovosti"
  ) +
  facet_wrap(~ vacc_status) +
  theme_minimal(base_size = 13)

4.1.4.2 Roční multipředpisovost

multi_presc_yearly <- multi_presc %>%
  mutate(year = year(date)) %>%
  filter(year > 2017)


multi_clients_yearly <- multi_presc_yearly %>%
  group_by(year, vacc_status, client_id) %>%
  summarise(presc_count = n(), .groups = "drop") %>%
  mutate(is_multi = presc_count > 1)


multi_summary_yearly <- multi_clients_yearly %>%
  group_by(year, vacc_status) %>%
  summarise(
    total_clients = n(),
    multi_clients = sum(is_multi),
    multi_rate = (multi_clients / total_clients) * 1000,
    .groups = "drop"
  )


multi_summary_yearly <- multi_summary_yearly %>%
  mutate(
    period_group = case_when(
      year < 2020 ~ "Pre-COVID",
      year > 2021 ~ "Post-COVID",
      TRUE ~ "COVID-Shock"
    )
  )


multi_summary_clean <- multi_summary_yearly %>%
  filter(period_group %in% c("Pre-COVID", "Post-COVID"))


t_test_results <- multi_summary_clean %>%
  group_by(vacc_status) %>%
  summarise(
    p_value = t.test(multi_rate ~ period_group)$p.value,
    mean_before = mean(multi_rate[period_group == "Pre-COVID"]),
    mean_after = mean(multi_rate[period_group == "Post-COVID"]),
    rate = mean_after / mean_before,
    .groups = "drop"
  )

V obou skupinách je možné sledovat nárůst. Ten však ani v jednom případě není statisticky významný. Zajímavé je však sledovat klesavou tendenci u neočkované populace mezi lety 2022 a 2023, zatímco u očkované populace je trend neustále rostoucí.

t_test_results %>%
  transmute(Vakcinace = vacc_status, p_value, `Průměr před` = mean_before, `Průměr po` = mean_after, `Poměr` = rate) %>%
  kable(format = "html", caption = "Porovnání roční multipředpisovosti v čase") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání roční multipředpisovosti v čase
Vakcinace p_value Průměr před Průměr po Poměr
Bez vakcinace 0.2265163 145.2588 193.5101 1.332175
Vakcinace 0.0917678 139.7487 183.4569 1.312763
ggplot(multi_summary_yearly, aes(x = year, y = multi_rate, color = vacc_status)) +
  geom_line(size = 0.6, alpha = 0.4) +
  stat_smooth(method = "loess", se = FALSE, size = 1.2) +
  annotate("rect",
           xmin = 2020, xmax = 2021,
           ymin = -Inf, ymax = Inf,
           alpha = 0.2, fill = "gray50") +
  scale_y_continuous(name = "Roční multipředpisovost na 1000 předepsaných klientů") +
  scale_x_continuous(breaks = unique(multi_summary_yearly$year)) +
  labs(
    title = "Roční multipředpisovost podle vakcinace",
    subtitle = "Šedá = covidové období (vyřazeno z t-testu)",
    x = "Rok",
    color = "Vakcinace"
  ) +
  theme_minimal()

Rozdíl v rozdílech nárůstů je také vysoko nad přijatelnou hladinou statistické významnosti.

did_data <- multi_summary_clean %>%
  mutate(
    vaccinated = if_else(vacc_status == "Vakcinace", 1, 0),
    post_covid = if_else(year >= 2022, 1, 0),
    interaction = vaccinated * post_covid
  )

did_model <- lm(log(multi_rate) ~ vaccinated * post_covid, data = did_data)
cat("```\n", paste(capture.output(summary(did_model)), collapse = "\n"), "\n```")
 
Call:
lm(formula = log(multi_rate) ~ vaccinated * post_covid, data = did_data)

Residuals:
       1        2        3        4        5        6        7        8 
-0.04517 -0.03659  0.04517  0.03659 -0.10415 -0.05379  0.10415  0.05379 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)            4.97750    0.06542  76.083 1.79e-07 ***
vaccinated            -0.03832    0.09252  -0.414   0.7000    
post_covid             0.28242    0.09252   3.053   0.0379 *  
vaccinated:post_covid -0.01106    0.13084  -0.085   0.9367    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.09252 on 4 degrees of freedom
Multiple R-squared:  0.8212,    Adjusted R-squared:  0.6871 
F-statistic: 6.123 on 3 and 4 DF,  p-value: 0.05625
 

4.1.4.3 Měsíční multipředpisovost

Neočkovaná skupina vykazuje po covidovém období téměř 34% nárůst. Růst je také znatelný ve skupině očkovaných, ten však není statisticky signifikantní.

t_test_results %>%
  transmute(Vakcinace = vacc_status, p_value, `Průměr před` = mean_before, `Průměr po` = mean_after, `Poměr` = rate) %>%
  kable(format = "html", caption = "Porovnání měsíční multipředpisovosti v čase") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání měsíční multipředpisovosti v čase
Vakcinace p_value Průměr před Průměr po Poměr
Bez vakcinace 0.0066103 78.75462 105.67934 1.341881
Vakcinace 0.5186993 89.80643 95.56067 1.064074
ggplot(multi_summary_monthly, aes(x = year_month, y = multi_rate, color = vacc_status)) +
  geom_line(size = 0.5, alpha = 0.4) + 
  stat_smooth(method = "loess", se = FALSE, size = 1.2) +  
  annotate("rect",
           xmin = as.Date("2020-01-01"), xmax = as.Date("2021-12-31"),
           ymin = -Inf, ymax = Inf, alpha = 0.2, fill = "gray50") +
  scale_y_continuous(name = "Měsíční multipředpisovost na 1000 předepsaných klientů") +
  labs(
    title = "Měsíční multipředpisovost podle vakcinačního statusu",
    subtitle = "Šedá = covidové období (vyřazeno z t-testu)",
    x = "Datum",
    color = "Vakcinace"
  ) +
  theme_minimal()

Interakce vakcíny a post-covidového období nevykazuje známky statistické významnosti.

did_data <- multi_summary_clean %>%
  mutate(
    vaccinated = if_else(vacc_status == "Vakcinace", 1, 0),
    post_covid = if_else(year_month >= as.Date("2022-01-01"), 1, 0),
    interaction = vaccinated * post_covid
  )

did_model <- lm(log(multi_rate) ~ vaccinated * post_covid, data = did_data)
cat("```\n", paste(capture.output(summary(did_model)), collapse = "\n"), "\n```")
 
Call:
lm(formula = log(multi_rate) ~ vaccinated * post_covid, data = did_data)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.94019 -0.19407  0.03713  0.26729  0.57878 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)            4.31147    0.07575  56.920   <2e-16 ***
vaccinated             0.12265    0.10712   1.145   0.2552    
post_covid             0.27837    0.10712   2.599   0.0109 *  
vaccinated:post_covid -0.20807    0.15149  -1.373   0.1729    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.3711 on 92 degrees of freedom
Multiple R-squared:  0.073, Adjusted R-squared:  0.04277 
F-statistic: 2.415 on 3 and 92 DF,  p-value: 0.07156
 

4.1.5 Způsob podání

Na první pohled je patrný trend zvyšování poměru injekcí a infuzí na úkor tablet. Injekce také mají všeobecně vyšší míru multipředpisovosti. To může částečně vysvětlovat všeobecně vyšší multipředpisovost v post-covidovém období. Otázkou je, zdali k nárůstu injekcí / infuzí dochází u dětí a mladistvých (s rostoucím věkem) přirozeně a je tedy možné vysvětlit změnu prostým stárnutím populace a nebo se jedná o dlouhodobu změnu.

prescriptions <- prescriptions %>%
  mutate(
    drug_form_std = case_when(
      is.na(drug_form) ~ NA_character_,
      str_detect(drug_form, "tableta|Tableta") ~ "Tableta",
      str_detect(drug_form, "tobolka|Tobolka") ~ "Tobolka",
      str_detect(drug_form, "injekční|Injekční|infuzní|Infuzní") ~ "Injekce / Infuze",
      str_detect(drug_form, "č[íi]pek|Čípek") ~ "Čípek",
      str_detect(drug_form, "perorální roztok|Perorální roztok") ~ "Perorální roztok",
      TRUE ~ "Other"
    )
  )

prescriptions <- prescriptions %>%
  mutate(
    year = year(presc_date)
  )

ggplot(prescriptions %>% filter(!is.na(presc_date), !is.na(drug_form)), aes(x = factor(year), fill = drug_form_std)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(
    title = "Poměr způsobu podání v čase",
    x = "Rok",
    y = "%",
    fill = "Způsob podání"
  ) +
  theme_minimal()

client_yearly_counts <- prescriptions %>%
  filter(!is.na(presc_date), !is.na(drug_form)) %>%
  group_by(client_id, year, drug_form_std) %>%
  summarise(presc_count = n(), .groups = "drop")

average_freq_per_year <- client_yearly_counts %>%
  group_by(year, drug_form_std) %>%
  summarise(avg_presc_per_client = mean(presc_count), .groups = "drop")

ggplot(average_freq_per_year, aes(x = factor(year), y = avg_presc_per_client, fill = drug_form_std)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(
    title = "Počet předpisů za rok na klienta dle způsobu podání",
    x = "Rok",
    y = "Průměrný počet předpisů na klienta",
    fill = "Způsob podání"
  ) +
  theme_minimal()

Poměr injekcí / infuzí na tablety v zásadě potvrzuje trend zobrazený v předchozím grafu.

filtered_data <- prescriptions %>%
  mutate(
    drug_form_std = ifelse(is.na(drug_form_std), "Unknown", drug_form_std)
  ) %>%
  filter(
    drug_form_std %in% c("Injekce / Infuze", "Tableta"),
    !is.na(vaccinated),
    !is.na(presc_yearmonth)
  )

form_vax_counts <- filtered_data %>%
  group_by(presc_yearmonth, vaccinated, drug_form_std) %>%
  summarise(n = n(), .groups = "drop")

month_totals <- filtered_data %>%
  group_by(presc_yearmonth, vaccinated) %>%
  summarise(total = n(), .groups = "drop")

form_vax_ratios <- form_vax_counts %>%
  left_join(month_totals, by = c("presc_yearmonth", "vaccinated")) %>%
  mutate(ratio = n / total)

form_vax_ratios <- form_vax_ratios %>%
  mutate(group = paste(drug_form_std, ifelse(vaccinated == "Vakcinace", "Vakcinace", "Bez vakcinace"), sep = " - "))

ggplot(form_vax_ratios, aes(x = presc_yearmonth, y = ratio, color = group)) +
  geom_line(size = 0.6) +                   
  geom_smooth(se = FALSE, span = 0.3) +   
  scale_y_continuous(labels = percent_format(accuracy = 1)) +
  scale_x_date(date_labels = "%Y",         
               date_breaks = "1 year",
               expand = c(0.01, 0)) +
  labs(
    title = "Poměr počtu předpisů podle způsobu podání a vakcinace",
    x = "Datum",
    y = "Poměr počtu předpisů",
    color = "Forma podání + Vakcinace"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0, hjust = 0.5))

4.2 Krátkodobý efekt

4.2.1 Prvopředpisy

4.2.1.1 SCCS vakcinovaných

Bohužel je patrné, že pro relevantní zhodnocení rizika prvopředpisu je v očkované populaci 5-11 příliš málo záznamů. Jediný významný výsledek je v pátem období (tedy 270 - 365 dní od vakcinace). Podle tohoto výsledku je riziko více než dvojnásobné (ovšem na velmi širokém intervalu spolehlivosti).

combined_vaccinated_prescribed <- dbGetQuery(con,
                                    "
                                    SELECT
                                    v.client_id,
                                    v.date AS vaccination_date,
                                    p.date AS event_date,
                                    'CPZP' as ins_company
                                    FROM cpzp_vaccinations v
                                    LEFT JOIN cpzp_clients c ON v.client_id = c.client_id
                                    LEFT JOIN cpzp_ins_start_end s ON c.client_id = s.client_id
                                    INNER JOIN cpzp_prescriptions p ON c.client_id = p.client_id
                                                                    AND p.presc_order = 1
                                                                    AND p.date BETWEEN v.date - INTERVAL '1 year' AND v.date + INTERVAL '1 year'
                                    WHERE v.event_order = 1
                                    AND s.ins_start_date < v.date - INTERVAL '1 year'
                                    AND s.ins_end_date > v.date + INTERVAL '1 year'
                                    AND v.date BETWEEN DATE '2021-08-04' - INTERVAL '7 days'
                                                       AND DATE '2021-08-04' + INTERVAL '7 days'
                                    AND 2021 - c.birthyear BETWEEN 12 AND 15
                                    
                                    UNION
                                    
                                    SELECT
                                    v.client_id,
                                    v.date AS vaccination_date,
                                    p.date AS event_date,
                                    'OZP' AS ins_company
                                    FROM OZP_vaccinations v
                                    LEFT JOIN ozp_clients c ON v.client_id = c.client_id
                                    LEFT JOIN ozp_ins_start_end s ON c.client_id = s.client_id
                                    INNER JOIN ozp_prescriptions p ON c.client_id = p.client_id
                                                                    AND p.presc_order = 1
                                                                    AND p.date BETWEEN v.date - INTERVAL '1 year' AND v.date + INTERVAL '1 year'
                                    WHERE v.event_order = 1
                                    AND s.ins_start_date < v.date - INTERVAL '1 year'
                                    AND s.ins_end_date > v.date + INTERVAL '1 year'
                                    AND v.date BETWEEN DATE '2021-08-04' - INTERVAL '7 days'
                                                       AND DATE '2021-08-04' + INTERVAL '7 days'
                                    AND 2021 - c.birthyear BETWEEN 12 AND 15
                                    ") %>%
  mutate(client_id = row_number()) %>%
  transmute(
    person_id = client_id,
    exposure_start_date = vaccination_date,
    event_date,
    observation_start_date = vaccination_date - days(365),
    observation_end_date = vaccination_date + days(365)
  )

sccs_vaccinated <- combined_vaccinated_prescribed %>%
  mutate(
    event_day = as.integer(event_date - observation_start_date),
    exposure_start_day = as.integer(exposure_start_date - observation_start_date),
    observation_start_day = 0,
    observation_end_day = as.integer(observation_end_date - observation_start_date)
  )

sccs_result_vaccinated <- standardsccs(
  formula = event_day ~ exposure_start_day,
  indiv = person_id,
  astart = observation_start_day,
  aend = observation_end_day,
  aevent = event_day,
  adrug = exposure_start_day,
  aedrug = exposure_start_day + 365,
  expogrp = list(c(0, 30, 90, 180, 270)),
  data = sccs_vaccinated
)

cat("```\n", paste(capture.output(print(sccs_result_vaccinated)), collapse = "\n"), "\n```")
 Call:
coxph(formula = Surv(rep(1, 480L), event) ~ exposure_start_day + 
    strata(indivL) + offset(log(interval)), data = chopdat, method = "exact")

  n= 480, number of events= 80 

                       coef exp(coef) se(coef)      z Pr(>|z|)   
exposure_start_day1  0.1008    1.1061   0.6030  0.167  0.86724   
exposure_start_day2 -0.3047    0.7374   0.5294 -0.575  0.56499   
exposure_start_day3  0.1008    1.1061   0.3761  0.268  0.78865   
exposure_start_day4  0.3885    1.4747   0.3371  1.152  0.24914   
exposure_start_day5  0.7835    2.1891   0.2880  2.721  0.00652 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                    exp(coef) exp(-coef) lower .95 upper .95
exposure_start_day1    1.1061     0.9041    0.3392     3.606
exposure_start_day2    0.7374     1.3562    0.2612     2.081
exposure_start_day3    1.1061     0.9041    0.5293     2.311
exposure_start_day4    1.4747     0.6781    0.7617     2.855
exposure_start_day5    2.1891     0.4568    1.2449     3.849

Concordance= 0.748  (se = 0.038 )
Likelihood ratio test= 8.52  on 5 df,   p=0.1
Wald test            = 9.2  on 5 df,   p=0.1
Score (logrank) test = 9.64  on 5 df,   p=0.09
 

4.2.1.2 SCCS nevakcinovaných

Neočkovaná část populace vykazuje mnohem relevantnější výsledky (i z důvodu většího množství sledovaných osob).

combined_unvaccinated_prescribed <- dbGetQuery(con,
                                               "
                                    SELECT
                                    p.client_id,
                                    DATE '2021-08-04' AS vaccination_date,
                                    p.date AS event_date,
                                    'CPZP' as ins_company
                                    FROM cpzp_prescriptions p
                                    LEFT JOIN cpzp_clients c ON P.client_id = c.client_id 
                                    LEFT JOIN cpzp_ins_start_end s ON c.client_id = s.client_id
                                    WHERE 1 = 1
                                    AND s.ins_start_date < DATE '2021-08-04'- INTERVAL '1 year'
                                    AND s.ins_end_date > DATE '2021-08-04' + INTERVAL '1 year'
                                    AND 2021 - c.birthyear BETWEEN 12 AND 15
                                    AND p.date BETWEEN DATE '2021-08-04' - INTERVAL '1 year' AND DATE '2021-08-04' + INTERVAL '1 year'
                                    AND p.presc_order = 1
                                    AND c.n_vacc = 0
                                    
                                    UNION
                                    
                                    SELECT
                                    p.client_id,
                                    DATE '2021-08-04' AS vaccination_date,
                                    p.date AS event_date,
                                    'OZP' as ins_company
                                    FROM ozp_prescriptions p
                                    LEFT JOIN ozp_clients c ON P.client_id = c.client_id 
                                    LEFT JOIN ozp_ins_start_end s ON c.client_id = s.client_id
                                    WHERE 1 = 1
                                    AND s.ins_start_date < DATE '2021-08-04' - INTERVAL '1 year'
                                    AND s.ins_end_date > DATE '2021-08-04' + INTERVAL '1 year'
                                    AND 2021 - c.birthyear BETWEEN 12 AND 15
                                    AND p.date BETWEEN DATE '2021-08-04' - INTERVAL '1 year' AND DATE '2021-08-04' + INTERVAL '1 year'
                                    AND p.presc_order = 1
                                    AND c.n_vacc = 0

                                               ") %>%
  mutate(client_id = row_number()) %>%
  transmute(
    person_id = client_id,
    exposure_start_date = vaccination_date,
    event_date,
    observation_start_date = vaccination_date - days(365),
    observation_end_date = vaccination_date + days(365)
  )

sccs_unvaccinated <- combined_unvaccinated_prescribed %>%
  mutate(
    event_day = as.integer(event_date - observation_start_date),
    exposure_start_day = as.integer(exposure_start_date - observation_start_date),
    observation_start_day = 0,
    observation_end_day = as.integer(observation_end_date - observation_start_date)
  )

sccs_result_unvaccinated <- standardsccs(
  formula = event_day ~ exposure_start_day,
  indiv = person_id,
  astart = observation_start_day,
  aend = observation_end_day,
  aevent = event_day,
  adrug = exposure_start_day,
  aedrug = exposure_start_day + 365,
  expogrp = list(c(0, 30, 90, 180, 270)),
  data = sccs_unvaccinated
)

cat("```\n", paste(capture.output(print(sccs_result_unvaccinated)), collapse = "\n"), "\n```")
 Call:
coxph(formula = Surv(rep(1, 4566L), event) ~ exposure_start_day + 
    strata(indivL) + offset(log(interval)), data = chopdat, method = "exact")

  n= 4566, number of events= 761 

                       coef exp(coef) se(coef)     z Pr(>|z|)   
exposure_start_day1 0.08217   1.08564  0.19380 0.424  0.67157   
exposure_start_day2 0.36553   1.44128  0.12674 2.884  0.00393 **
exposure_start_day3 0.27970   1.32274  0.11185 2.501  0.01240 * 
exposure_start_day4 0.32579   1.38513  0.10994 2.963  0.00304 **
exposure_start_day5 0.27911   1.32196  0.10921 2.556  0.01060 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                    exp(coef) exp(-coef) lower .95 upper .95
exposure_start_day1     1.086     0.9211    0.7425     1.587
exposure_start_day2     1.441     0.6938    1.1243     1.848
exposure_start_day3     1.323     0.7560    1.0623     1.647
exposure_start_day4     1.385     0.7220    1.1166     1.718
exposure_start_day5     1.322     0.7565    1.0672     1.637

Concordance= 0.709  (se = 0.014 )
Likelihood ratio test= 17.88  on 5 df,   p=0.003
Wald test            = 17.87  on 5 df,   p=0.003
Score (logrank) test = 18.02  on 5 df,   p=0.003
 

4.2.1.3 Porovnání rizik

Podle naměřených hodnot je riziko u očkovaných osob mírně nižší nebo srovnatelné až 180 dní od vakcinace, poté se naopak mírně zvýší. P-hodnoty rozdílů mezi skupinami jsou však ve všech případech vysoko nad akceptovatelnou úrovní (což je i mimo výpočet patrné z širokých intervalů spolehlivosti).

coef_unvax <- sccs_result_unvaccinated$coefficients
conf_unvax <- sccs_result_unvaccinated$conf.int

df_unvax <- data.frame(
  term = rownames(coef_unvax),
  logHR = coef_unvax[, "coef"],
  SE = coef_unvax[, "se(coef)"],
  HR = coef_unvax[, "exp(coef)"],
  conf.low = conf_unvax[, "lower .95"],
  conf.high = conf_unvax[, "upper .95"]
)

coef_vax <- sccs_result_vaccinated$coefficients
conf_vax <- sccs_result_vaccinated$conf.int

df_vax <- data.frame(
  term = rownames(coef_vax),
  logHR = coef_vax[, "coef"],
  SE = coef_vax[, "se(coef)"],
  HR = coef_vax[, "exp(coef)"],
  conf.low = conf_vax[, "lower .95"],
  conf.high = conf_vax[, "upper .95"]
)

df_unvax <- df_unvax %>% mutate(group = "Bez Vakcinace")
df_vax <- df_vax %>% mutate(group = "Vakcinace")

combined <- bind_rows(df_vax, df_unvax)

combined$period <- as.integer(gsub("exposure_start_day", "", combined$term))

start_days <- c(0, 30, 90, 180, 270)
end_days   <- c(29, 89, 179, 269, Inf)

period_labels <- ifelse(
  is.infinite(end_days),
  paste0(start_days, "+"),
  paste0(start_days, "-", end_days)
)

combined$period_label <- factor(combined$period,
                                levels = seq_along(period_labels),
                                labels = period_labels)

ggplot(combined, aes(x = period_label, y = HR, color = group)) +
  geom_point(position = position_dodge(width = 0.6), size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2,
                position = position_dodge(width = 0.6)) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  labs(
    x = "Období expozice (dny)",
    y = "Poměr rizik (HR)",
    title = "Porovnání SCCS: očkovaní vs neočkovaní"
  ) +
  theme_minimal()

vax_df <- combined %>%
  filter(group == "Vakcinace") %>%
  rename_with(~ paste0(., "_vax"), -c(period, period_label))

unvax_df <- combined %>%
  filter(group == "Bez Vakcinace") %>%
  rename_with(~ paste0(., "_unvax"), -c(period, period_label))

df_diff <- inner_join(vax_df, unvax_df, by = c("period", "period_label"))

df_diff <- df_diff %>%
  mutate(
    logHR_diff = logHR_vax - logHR_unvax,
    se_diff = sqrt(SE_vax^2 + SE_unvax^2),
    ci_low = logHR_diff - 1.96 * se_diff,
    ci_high = logHR_diff + 1.96 * se_diff
  )

ggplot(df_diff, aes(x = period_label, y = logHR_diff)) +
  geom_point(color = "blue", size = 3) +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.2, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "Rozdíl v log(HR) (Očkovaní - Neočkovaní)",
    x = "Období expozice",
    y = "Rozdíl v log(HR)"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

common_terms <- intersect(df_vax$term, df_unvax$term)

did_results <- data.frame(
  term = character(),
  period_label = character(),
  logHR_vax = numeric(),
  logHR_unvax = numeric(),
  DiD_logHR = numeric(),
  DiD_SE = numeric(),
  Z = numeric(),
  P_value = numeric(),
  stringsAsFactors = FALSE
)

for (term in common_terms) {
  logHR_vax <- df_vax %>% filter(term == !!term) %>% pull(logHR)
  SE_vax    <- df_vax %>% filter(term == !!term) %>% pull(SE)
  
  logHR_unvax <- df_unvax %>% filter(term == !!term) %>% pull(logHR)
  SE_unvax    <- df_unvax %>% filter(term == !!term) %>% pull(SE)
  
  DiD_logHR <- logHR_vax - logHR_unvax
  DiD_SE <- sqrt(SE_vax^2 + SE_unvax^2)
  Z <- DiD_logHR / DiD_SE
  P_value <- 2 * pnorm(-abs(Z))
  
  period_idx <- as.integer(gsub("exposure_start_day", "", term))
  period_label <- period_labels[period_idx]
  
  did_results <- rbind(did_results, data.frame(
    term = term,
    period_label = period_label,
    logHR_vax = logHR_vax,
    logHR_unvax = logHR_unvax,
    DiD_logHR = DiD_logHR,
    DiD_SE = DiD_SE,
    Z = Z,
    P_value = P_value
  ))
}

did_results %>%
  transmute(Období = period_label, `log(HR) Vakcinace` = logHR_vax, `log(HR) Bez vakcinace` = logHR_unvax, `DiD log(HR)` = DiD_logHR, p_value = P_value) %>%
  kable(format = "html", caption = "Porovnání rizik očkovaných a neočkovaných (SCCS)") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání rizik očkovaných a neočkovaných (SCCS)
Období log(HR) Vakcinace log(HR) Bez vakcinace DiD log(HR) p_value
0-29 0.1008047 0.0821706 0.0186341 0.9765303
30-89 -0.3046604 0.3655330 -0.6701934 0.2182940
90-179 0.1008047 0.2797016 -0.1788969 0.6484026
180-269 0.3884868 0.3257927 0.0626941 0.8596533
270+ 0.7834806 0.2791118 0.5043688 0.1015081

4.2.2 Předpisy celkem

4.2.2.1 SCCS vakcinovaných

U celkového počtu předpisů je výhoda ve větším množství záznamů a tedy i vyšší pravděpodobnosti, že některý z výsledků bude významný. To se bohužel u očkované skupiny nepotvrdilo a jediný signifikantní výsledek je opět v období 270 - 365 dní po vakcinaci.

combined_vaccinated_prescribed <- dbGetQuery(con,
                                    "
                                    SELECT
                                    v.client_id,
                                    v.date AS vaccination_date,
                                    p.date AS event_date,
                                    d.drug_form,
                                    'CPZP' as ins_company
                                    FROM cpzp_vaccinations v
                                    LEFT JOIN cpzp_clients c ON v.client_id = c.client_id
                                    LEFT JOIN cpzp_ins_start_end s ON c.client_id = s.client_id
                                    INNER JOIN cpzp_prescriptions p ON c.client_id = p.client_id
                                                                    AND p.date BETWEEN v.date - INTERVAL '1 year' AND v.date + INTERVAL '1 year'
                                    LEFT JOIN cpzp_drugs_catalog d ON p.presc_id  = d.presc_id 
                                    WHERE v.event_order = 1
                                    AND s.ins_start_date < v.date - INTERVAL '1 year'
                                    AND s.ins_end_date > v.date + INTERVAL '1 year'
                                    AND v.date BETWEEN DATE '2021-08-04' - INTERVAL '7 days'
                                                       AND DATE '2021-08-04' + INTERVAL '7 days'
                                    AND 2021 - c.birthyear BETWEEN 12 AND 15
                                    
                                    UNION
                                    
                                    SELECT
                                    v.client_id,
                                    v.date AS vaccination_date,
                                    p.date AS event_date,
                                    d.drug_form,
                                    'OZP' AS ins_company
                                    FROM OZP_vaccinations v
                                    LEFT JOIN ozp_clients c ON v.client_id = c.client_id
                                    LEFT JOIN ozp_ins_start_end s ON c.client_id = s.client_id
                                    INNER JOIN ozp_prescriptions p ON c.client_id = p.client_id
                                                                    AND p.date BETWEEN v.date - INTERVAL '1 year' AND v.date + INTERVAL '1 year'
                                    LEFT JOIN ozp_drugs_catalog d ON p.presc_id  = d.presc_id 
                                    WHERE v.event_order = 1
                                    AND s.ins_start_date < v.date - INTERVAL '1 year'
                                    AND s.ins_end_date > v.date + INTERVAL '1 year'
                                    AND v.date BETWEEN DATE '2021-08-04' - INTERVAL '7 days'
                                                       AND DATE '2021-08-04' + INTERVAL '7 days'
                                    AND 2021 - c.birthyear BETWEEN 12 AND 15
                                    ") %>%
  mutate(client_id = as.numeric(factor(client_id)),
         drug_form = tolower(drug_form),
         inj_or_inf = str_detect(drug_form, "injekční|infuzní"))

injectable_infusional <- combined_vaccinated_prescribed %>%
  filter(inj_or_inf) %>%
  arrange(client_id, event_date) %>%
  group_by(client_id) %>%
  mutate(
    day_diff = as.numeric(difftime(event_date, lag(event_date, default = first(event_date)), units = "days")),
    new_window = is.na(day_diff) | day_diff > 10,
    window_id = cumsum(new_window)
  ) %>%
  group_by(client_id, window_id) %>%
  slice(1) %>%
  ungroup()

other_drugs <- combined_vaccinated_prescribed %>%
  filter(!inj_or_inf)

combined_vaccinated_prescribed <- bind_rows(injectable_infusional, other_drugs) %>%
  arrange(client_id, event_date) %>%
  transmute(
    person_id = client_id,
    exposure_start_date = vaccination_date,
    event_date,
    observation_start_date = vaccination_date - days(365),
    observation_end_date = vaccination_date + days(365)
  )

sccs_vaccinated <- combined_vaccinated_prescribed %>%
  mutate(
    event_day = as.integer(event_date - observation_start_date),
    exposure_start_day = as.integer(exposure_start_date - observation_start_date),
    observation_start_day = 0,
    observation_end_day = as.integer(observation_end_date - observation_start_date)
  )

sccs_result_vaccinated <- standardsccs(
  formula = event_day ~ exposure_start_day,
  indiv = person_id,
  astart = observation_start_day,
  aend = observation_end_day,
  aevent = event_day,
  adrug = exposure_start_day,
  aedrug = exposure_start_day + 365,
  expogrp = list(c(0, 30, 90, 180, 270)),
  data = sccs_vaccinated
)

cat("```\n", paste(capture.output(print(sccs_result_vaccinated)), collapse = "\n"), "\n```")
 Call:
coxph(formula = Surv(rep(1, 1122L), event) ~ exposure_start_day + 
    strata(indivL) + offset(log(interval)), data = chopdat, method = "exact")

  n= 1122, number of events= 187 

                        coef exp(coef) se(coef)      z Pr(>|z|)   
exposure_start_day1 -0.14036   0.86905  0.42258 -0.332  0.73978   
exposure_start_day2  0.01379   1.01389  0.28868  0.048  0.96189   
exposure_start_day3 -0.08629   0.91733  0.25404 -0.340  0.73410   
exposure_start_day4  0.18815   1.20701  0.22783  0.826  0.40890   
exposure_start_day5  0.56829   1.76525  0.19377  2.933  0.00336 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                    exp(coef) exp(-coef) lower .95 upper .95
exposure_start_day1    0.8690     1.1507    0.3796     1.989
exposure_start_day2    1.0139     0.9863    0.5758     1.785
exposure_start_day3    0.9173     1.0901    0.5576     1.509
exposure_start_day4    1.2070     0.8285    0.7723     1.886
exposure_start_day5    1.7653     0.5665    1.2075     2.581

Concordance= 0.752  (se = 0.026 )
Likelihood ratio test= 9.6  on 5 df,   p=0.09
Wald test            = 10.51  on 5 df,   p=0.06
Score (logrank) test = 10.79  on 5 df,   p=0.06
 

4.2.2.2 SCCS nevakcinovaných

Naopak u neočkované populace je díky výrazně většímu množství záznamů možné zhodnotit rizika s větší jistotou. Výsledek v zásadě říká, že po celou dobu od rozhodného data (vyjma prvního měsíce) je riziko zhruba jeden a půl násobné.

combined_unvaccinated_prescribed <- dbGetQuery(con,
                                               "
                                    SELECT
                                    p.client_id,
                                    DATE '2021-08-04' AS vaccination_date,
                                    p.date AS event_date,
                                    d.drug_form,
                                    'CPZP' as ins_company
                                    FROM cpzp_prescriptions p
                                    LEFT JOIN cpzp_clients c ON P.client_id = c.client_id 
                                    LEFT JOIN cpzp_ins_start_end s ON c.client_id = s.client_id
                                    LEFT JOIN cpzp_drugs_catalog d ON p.presc_id  = d.presc_id 
                                    WHERE 1 = 1
                                    AND s.ins_start_date < DATE '2021-08-04' - INTERVAL '1 year'
                                    AND s.ins_end_date > DATE '2021-08-04' + INTERVAL '1 year'
                                    AND 2021 - c.birthyear BETWEEN 12 AND 15
                                    AND p.date BETWEEN DATE '2021-08-04' - INTERVAL '1 year' AND DATE '2021-08-04' + INTERVAL '1 year'
                                    AND c.n_vacc = 0
                                    
                                    UNION
                                    
                                    SELECT
                                    p.client_id,
                                    DATE '2021-08-04' AS vaccination_date,
                                    p.date AS event_date,
                                    d.drug_form,
                                    'OZP' as ins_company
                                    FROM ozp_prescriptions p
                                    LEFT JOIN ozp_clients c ON P.client_id = c.client_id 
                                    LEFT JOIN ozp_ins_start_end s ON c.client_id = s.client_id
                                    LEFT JOIN ozp_drugs_catalog d ON p.presc_id  = d.presc_id 
                                    WHERE 1 = 1
                                    AND s.ins_start_date < DATE '2021-08-04' - INTERVAL '1 year'
                                    AND s.ins_end_date > DATE '2021-08-04' + INTERVAL '1 year'
                                    AND 2021 - c.birthyear BETWEEN 12 AND 15
                                    AND p.date BETWEEN DATE '2021-08-04' - INTERVAL '1 year' AND DATE '2021-08-04' + INTERVAL '1 year'
                                    AND c.n_vacc = 0

                                               ") %>%
  mutate(client_id = as.numeric(factor(client_id)),
         drug_form = tolower(drug_form),
         inj_or_inf = str_detect(drug_form, "injekční|infuzní"))

injectable_infusional <- combined_unvaccinated_prescribed %>%
  filter(inj_or_inf) %>%
  arrange(client_id, event_date) %>%
  group_by(client_id) %>%
  mutate(
    day_diff = as.numeric(difftime(event_date, lag(event_date, default = first(event_date)), units = "days")),
    new_window = is.na(day_diff) | day_diff > 10,
    window_id = cumsum(new_window)
  ) %>%
  group_by(client_id, window_id) %>%
  slice(1) %>%
  ungroup()

other_drugs <- combined_unvaccinated_prescribed %>%
  filter(!inj_or_inf)

combined_unvaccinated_prescribed <- bind_rows(injectable_infusional, other_drugs) %>%
  arrange(client_id, event_date) %>%
  transmute(
    person_id = client_id,
    exposure_start_date = vaccination_date,
    event_date,
    observation_start_date = vaccination_date - days(365),
    observation_end_date = vaccination_date + days(365)
  )

sccs_unvaccinated <- combined_unvaccinated_prescribed %>%
  mutate(
    event_day = as.integer(event_date - observation_start_date),
    exposure_start_day = as.integer(exposure_start_date - observation_start_date),
    observation_start_day = 0,
    observation_end_day = as.integer(observation_end_date - observation_start_date)
  )

sccs_result_unvaccinated <- standardsccs(
  formula = event_day ~ exposure_start_day,
  indiv = person_id,
  astart = observation_start_day,
  aend = observation_end_day,
  aevent = event_day,
  adrug = exposure_start_day,
  aedrug = exposure_start_day + 365,
  expogrp = list(c(0, 30, 90, 180, 270)),
  data = sccs_unvaccinated
)

cat("```\n", paste(capture.output(print(sccs_result_unvaccinated)), collapse = "\n"), "\n```")
 Call:
coxph(formula = Surv(rep(1, 10308L), event) ~ exposure_start_day + 
    strata(indivL) + offset(log(interval)), data = chopdat, method = "exact")

  n= 10308, number of events= 1718 

                       coef exp(coef) se(coef)     z Pr(>|z|)    
exposure_start_day1 0.08648   1.09033  0.13151 0.658    0.511    
exposure_start_day2 0.39769   1.48838  0.08507 4.675 2.94e-06 ***
exposure_start_day3 0.33367   1.39608  0.07453 4.477 7.57e-06 ***
exposure_start_day4 0.43198   1.54030  0.07189 6.009 1.87e-09 ***
exposure_start_day5 0.38231   1.46567  0.07150 5.347 8.95e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                    exp(coef) exp(-coef) lower .95 upper .95
exposure_start_day1     1.090     0.9172    0.8426     1.411
exposure_start_day2     1.488     0.6719    1.2598     1.758
exposure_start_day3     1.396     0.7163    1.2063     1.616
exposure_start_day4     1.540     0.6492    1.3379     1.773
exposure_start_day5     1.466     0.6823    1.2740     1.686

Concordance= 0.705  (se = 0.009 )
Likelihood ratio test= 63.17  on 5 df,   p=3e-12
Wald test            = 62.83  on 5 df,   p=3e-12
Score (logrank) test = 63.63  on 5 df,   p=2e-12
 

4.2.2.3 Porovnání rizik

Rozdíl mezi riziky je všeobecně vyhodnocen jako menší pro očkovanou populaci. Ani jedno z porovnání však není statisticky významné.

coef_unvax <- sccs_result_unvaccinated$coefficients
conf_unvax <- sccs_result_unvaccinated$conf.int

df_unvax <- data.frame(
  term = rownames(coef_unvax),
  logHR = coef_unvax[, "coef"],
  SE = coef_unvax[, "se(coef)"],
  HR = coef_unvax[, "exp(coef)"],
  conf.low = conf_unvax[, "lower .95"],
  conf.high = conf_unvax[, "upper .95"]
)

coef_vax <- sccs_result_vaccinated$coefficients
conf_vax <- sccs_result_vaccinated$conf.int

df_vax <- data.frame(
  term = rownames(coef_vax),
  logHR = coef_vax[, "coef"],
  SE = coef_vax[, "se(coef)"],
  HR = coef_vax[, "exp(coef)"],
  conf.low = conf_vax[, "lower .95"],
  conf.high = conf_vax[, "upper .95"]
)

df_unvax <- df_unvax %>% mutate(group = "Unvaccinated")
df_vax <- df_vax %>% mutate(group = "Vaccinated")

combined <- bind_rows(df_vax, df_unvax)

combined$period <- as.integer(gsub("exposure_start_day", "", combined$term))

start_days <- c(0, 30, 90, 180, 270)
end_days   <- c(29, 89, 179, 269, Inf)

period_labels <- ifelse(
  is.infinite(end_days),
  paste0(start_days, "+"),
  paste0(start_days, "-", end_days)
)

combined$period_label <- factor(combined$period,
                                levels = seq_along(period_labels),
                                labels = period_labels)

ggplot(combined, aes(x = period_label, y = HR, color = group)) +
  geom_point(position = position_dodge(width = 0.6), size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2,
                position = position_dodge(width = 0.6)) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  labs(
    x = "Období expozice (dny)",
    y = "Poměr rizik (HR)",
    title = "Porovnání SCCS: očkovaní vs neočkovaní"
  ) +
  theme_minimal()

vax_df <- combined %>%
  filter(group == "Vaccinated") %>%
  rename_with(~ paste0(., "_vax"), -c(period, period_label))

unvax_df <- combined %>%
  filter(group == "Unvaccinated") %>%
  rename_with(~ paste0(., "_unvax"), -c(period, period_label))

df_diff <- inner_join(vax_df, unvax_df, by = c("period", "period_label"))

df_diff <- df_diff %>%
  mutate(
    logHR_diff = logHR_vax - logHR_unvax,
    se_diff = sqrt(SE_vax^2 + SE_unvax^2),
    ci_low = logHR_diff - 1.96 * se_diff,
    ci_high = logHR_diff + 1.96 * se_diff
  )

ggplot(df_diff, aes(x = period_label, y = logHR_diff)) +
  geom_point(color = "blue", size = 3) +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.2, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "Rozdíl v log(HR) (Očkovaní - Neočkovaní)",
    x = "Období expozice",
    y = "Rozdíl v log(HR)"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

common_terms <- intersect(df_vax$term, df_unvax$term)

did_results <- data.frame(
  term = character(),
  period_label = character(),
  logHR_vax = numeric(),
  logHR_unvax = numeric(),
  DiD_logHR = numeric(),
  DiD_SE = numeric(),
  Z = numeric(),
  P_value = numeric(),
  stringsAsFactors = FALSE
)

for (term in common_terms) {
  logHR_vax <- df_vax %>% filter(term == !!term) %>% pull(logHR)
  SE_vax    <- df_vax %>% filter(term == !!term) %>% pull(SE)
  
  logHR_unvax <- df_unvax %>% filter(term == !!term) %>% pull(logHR)
  SE_unvax    <- df_unvax %>% filter(term == !!term) %>% pull(SE)
  
  DiD_logHR <- logHR_vax - logHR_unvax
  DiD_SE <- sqrt(SE_vax^2 + SE_unvax^2)
  Z <- DiD_logHR / DiD_SE
  P_value <- 2 * pnorm(-abs(Z))
  
  period_idx <- as.integer(gsub("exposure_start_day", "", term))
  period_label <- period_labels[period_idx]
  
  did_results <- rbind(did_results, data.frame(
    term = term,
    period_label = period_label,
    logHR_vax = logHR_vax,
    logHR_unvax = logHR_unvax,
    DiD_logHR = DiD_logHR,
    DiD_SE = DiD_SE,
    Z = Z,
    P_value = P_value
  ))
}

did_results %>%
  transmute(Období = period_label, `log(HR) Vakcinace` = logHR_vax, `log(HR) Bez vakcinace` = logHR_unvax, `DiD log(HR)` = DiD_logHR, p_value = P_value) %>%
  kable(format = "html", caption = "Porovnání rizik (SCCS) očkovaných a neočkovaných") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání rizik (SCCS) očkovaných a neočkovaných
Období log(HR) Vakcinace log(HR) Bez vakcinace DiD log(HR) p_value
0-29 -0.1403574 0.0864778 -0.2268352 0.6082714
30-89 0.0137933 0.3976904 -0.3838971 0.2020880
90-179 -0.0862901 0.3336685 -0.4199587 0.1126796
180-269 0.1881467 0.4319795 -0.2438327 0.3074186
270+ 0.5682940 0.3823111 0.1859829 0.3678679

4.3 Shrnutí

Z dlouhodobého hlediska je nárůst ve všech sledovaných kategoriích naprosto zřejmý. Rozdíl mezi očkovanou a neočkovanou populací však není možné prokázat ani v jednom případě. Změny ve způsobu podání (vyšší podíl injekcí a infuzí na úkor tablet) mohou částečně vysvětlit zvýšenou multipředpisovost, je však otázkou zdali se jedná o přirozený jev způsobený stárnutím populace a nebo je vliv jiný.

Co se týče vyhodnocení rizika v krátkodobém horizontu, trpí dataset na nedostatek relevantních dat. Ve sledovaných obdobích nebyl nalezen jediný statisticky významný rozdíl mezi očkovanou a neočkovanou populací.

5 Věková kategorie 16-29

5.1 Dlouhodobý efekt

5.1.1 Distribuce věku a pohlaví

I v této věkové kategorii jsou patrné rozdíly, výrazný rozdíl způsobený odlišnými riziky na základě věku by však ani zde neměl nastat.

unique_clients <- prescriptions %>%
  filter(!is.na(age_2021), !is.na(sex), !is.na(vaccinated)) %>%
  distinct(client_id, .keep_all = TRUE)  

age_data <- unique_clients %>%
  group_by(vaccinated, age_2021) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(vaccinated) %>%
  mutate(prop = count / sum(count),
         variable = "Age",
         value = as.character(age_2021)) %>%
  ungroup()

sex_data <- unique_clients %>%
  group_by(vaccinated, sex) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(vaccinated) %>%
  mutate(prop = count / sum(count),
         variable = "Sex",
         value = sex) %>%
  ungroup()
ggplot(unique_clients, aes(x = age_2021, fill = factor(vaccinated))) +
  geom_histogram(aes(y = after_stat(density)), position = "identity", alpha = 0.4, binwidth = 1) +
  scale_fill_brewer(palette = "Set1", name = "Vaccinated") +
  scale_x_continuous(breaks = seq(floor(min(unique_clients$age)), ceiling(max(unique_clients$age)), by = 1)) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    title = "Věkové rozložení podle statusu vakcinace",
    x = "Věk", y = "%",
    fill = "Vakcinace"
  ) +
  theme_minimal()

ggplot(sex_data, aes(x = sex, y = prop, color = vaccinated, group = vaccinated)) +
  geom_line(size = 1) +
  geom_point(size = 2) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    title = "Pohlaví podle vakcinace",
    x = "Pohlaví",
    y = "%",
    color = "Vakcinace"
  ) +
  theme_minimal()

5.1.2 Předpověď vs realita

prescriptions_clean <- prescriptions %>%
  filter(!is.na(presc_date)) %>%
  mutate(presc_date = as.Date(presc_date)) %>%
  arrange(presc_date) %>%
  filter(presc_date >= as.Date("2015-01-01") & presc_date <= as.Date("2023-12-31")) %>%
  mutate(count = 1)  


results_vacc <- prescriptions_clean %>%
  filter(!is.na(vaccinated)) %>%
  group_by(vaccinated) %>%
  group_split() %>%
  map(~ {
    group_val <- unique(.x$vaccinated)
    analyze_group(df = filter(prescriptions_clean, vaccinated == group_val), 
                  group_id = group_val, 
                  group_type = "Vaccination")
  }) %>%
  compact() %>%
  bind_rows()

results_sex <- prescriptions_clean %>%
  filter(!is.na(sex)) %>%
  group_by(sex) %>%
  group_split() %>%
  map(~ {
    group_val <- unique(.x$sex)
    analyze_group(df = filter(prescriptions_clean, sex == group_val), 
                  group_id = group_val, 
                  group_type = "Sex")
  }) %>%
  compact() %>%
  bind_rows()

results_all <- bind_rows(results_vacc, results_sex) %>%
  arrange(p_value)

plot_data_all <- results_all %>%
  unnest(plot_data, names_sep = "_")

Výkyvy od předpokládaného trendu jsou zde ve všech skupinách. Vyjma vakcinované kategorie jsou také všechny statisticky významné.

results_all %>%
  mutate(Skupina = group, RMSE = rmse, rRMSE = rrmse) %>%
  select(
    Skupina, RMSE, rRMSE, t_stat, p_value
  ) %>%
  kable(format = "html", caption = "Porovnání předpisovosti v čase (ETS)") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání předpisovosti v čase (ETS)
Skupina RMSE rRMSE t_stat p_value
F 140.34671 0.2687347 12.2492343 0.0000000
M 115.36129 0.2724855 -6.5066174 0.0000000
Bez vakcinace 40.72322 0.1418763 2.7188140 0.0091492
Vakcinace 96.13567 0.1459734 -0.1575997 0.8754474
ggplot(plot_data_all, aes(x = plot_data_month)) +
  geom_line(aes(y = plot_data_prescriptions), color = "blue") +
  geom_line(aes(y = plot_data_predicted), color = "red", linetype = "dashed") +
  geom_vline(xintercept = as.Date("2020-01-01"), linetype = "dotted") +
  facet_wrap(~ paste(plot_data_group_type, plot_data_group, sep = ": "), scales = "free_y") +
  labs(
    title = "Realita vs Předpověď (ETS)",
    subtitle = "Modrá = Realita, Červená = Předpověď",
    x = "Datum", y = "Počet předpisů"
  ) +
  theme_minimal()

5.1.3 Prvopředpisovost

first_presc <- dbGetQuery(con,
                          "
                          SELECT
                            c.client_id,
                            c.n_vacc,
                            c.n_presc,
                            c.sex,
                            c.birthyear,
                            FIRST(p.date) OVER (PARTITION BY c.client_id ORDER BY p.date ASC) AS first_presc_date,
                            'CPZP' AS ins_company
                          FROM
                            cpzp_clients c
                          LEFT JOIN cpzp_ins_start_end i ON c.client_id = i.client_id
                          LEFT JOIN cpzp_prescriptions p ON c.client_id = p.client_id AND p.presc_order = 1
                          WHERE
                            i.ins_start_date < '2015-01-01' AND i.ins_end_date = '9999-12-31'
                            AND 2021 - c.birthyear BETWEEN 16 AND 29
                            
                          UNION
                          
                          SELECT
                            c.client_id,
                            c.n_vacc,
                            c.n_presc,
                            c.sex,
                            c.birthyear,
                            FIRST(p.date) OVER (PARTITION BY c.client_id ORDER BY p.date ASC) AS first_presc_date,
                            'OZP' AS ins_company
                          FROM
                            ozp_clients c
                          LEFT JOIN ozp_ins_start_end i ON c.client_id = i.client_id
                          LEFT JOIN ozp_prescriptions p ON c.client_id = p.client_id AND p.presc_order = 1
                          WHERE
                            i.ins_start_date < '2015-01-01' AND i.ins_end_date IS NULL
                            AND 2021 - c.birthyear BETWEEN 16 AND 29
                          ") %>%
  mutate(first_presc_yearmonth = floor_date(first_presc_date, unit = "month"),
         sex = ifelse(sex == "Z", "F", sex),
         vacc_status = ifelse(n_vacc > 0, "Vakcinace", "Bez vakcinace"),
         age_2021 = 2021 - birthyear)

monthly_first_presc <- first_presc %>%
  count(first_presc_yearmonth, name = "new_prescriptions") %>%
  arrange(first_presc_yearmonth) %>%
  mutate(
    cum_new_presc = cumsum(new_prescriptions),  
    total_population = nrow(first_presc),
    at_risk_population = total_population - lag(cum_new_presc, default = 0),
    rate = new_prescriptions / at_risk_population
  ) %>%
  mutate(
    period_group = case_when(
      first_presc_yearmonth < as.Date("2020-01-01") ~ "Pre-COVID",
      first_presc_yearmonth > as.Date("2021-12-31") ~ "Post-COVID",
      TRUE ~ "COVID-Shock"
    )
  )

monthly_clean <- monthly_first_presc %>%
  filter(first_presc_yearmonth >= as.Date("2018-01-01")) %>%
  filter(period_group %in% c("Pre-COVID", "Post-COVID"))

ttest <- t.test(rate ~ period_group, data = monthly_clean)

V porovnání před-covidového a po-covidového období došlo k jednoznačnému nárůstu prvopředpisovosti.

ggplot(monthly_first_presc %>% filter(first_presc_yearmonth >= as.Date("2018-01-01")),
       aes(x = first_presc_yearmonth, y = rate * 1000)) +
  geom_line(color = "steelblue", size = 1) +
  annotate("rect",
           xmin = as.Date("2020-01-01"), xmax = as.Date("2021-12-31"),
           ymin = -Inf, ymax = Inf, alpha = 0.2, fill = "gray50") +
  scale_y_continuous(name = "Počet prvopředpisů na 1000 osob bez předpisu", limits = c(0, 3)) +
  labs(
    title = "Prvopředpisovost v čase",
    subtitle = "Šedá = covidové období (vyřazeno z t-testu)",
    x = "Datum"
  ) +
  theme_minimal()

cat("```\n", paste(capture.output(ttest), collapse = "\n"), "\n```")
 
    Welch Two Sample t-test

data:  rate by period_group
t = 5.6304, df = 45.941, p-value = 1.037e-06
alternative hypothesis: true difference in means between group Post-COVID and group Pre-COVID is not equal to 0
95 percent confidence interval:
 0.0002250338 0.0004754821
sample estimates:
mean in group Post-COVID  mean in group Pre-COVID 
             0.002210366              0.001860108 
 
total_by_vacc <- first_presc %>%
  count(vacc_status, name = "total_population")

monthly_by_vacc <- first_presc %>%
  count(vacc_status, first_presc_yearmonth, name = "new_prescriptions") %>%
  arrange(vacc_status, first_presc_yearmonth) %>%
  group_by(vacc_status) %>%
  mutate(
    cum_prescriptions = cumsum(new_prescriptions),
    total_population = total_by_vacc$total_population[match(vacc_status, total_by_vacc$vacc_status)],
    at_risk_population = total_population - lag(cum_prescriptions, default = 0),
    rate = new_prescriptions / at_risk_population
  ) %>%
  ungroup() %>%
  mutate(
    period_group = case_when(
      first_presc_yearmonth < as.Date("2020-01-01") ~ "Pre-COVID",
      first_presc_yearmonth > as.Date("2021-12-31") ~ "Post-COVID",
      TRUE ~ "COVID-Shock"
    )
  )
  
monthly_by_vacc_clean <- monthly_by_vacc %>%
  filter(first_presc_yearmonth >= as.Date("2018-01-01")) %>%
  filter(period_group %in% c("Pre-COVID", "Post-COVID"))

first_presc_table <- monthly_by_vacc_clean %>%
  group_by(vacc_status) %>%
  summarise(
    p_value = t.test(rate ~ period_group)$p.value,
    mean_before = mean(rate[period_group == "Pre-COVID"]),
    mean_after = mean(rate[period_group == "Post-COVID"]),
    rate = mean_after / mean_before
  )

Rozdíly mezi skupinami jsou však nepatrné až téměř nulové. To také potvrzuje p-hodnota příslušného DiD modelu.

first_presc_table %>%
  transmute(Vakcinace = vacc_status, p_value, `Průměr před` = mean_before * 1000, `Průměr po` = mean_after * 1000, `Poměr` = rate) %>%
  kable(format = "html", caption = "Porovnání prvopředpisovosti v čase podle vakcinace") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání prvopředpisovosti v čase podle vakcinace
Vakcinace p_value Průměr před Průměr po Poměr
Bez vakcinace 0.0001376 1.707657 2.032872 1.190445
Vakcinace 0.0000011 1.934736 2.298202 1.187863
ggplot(monthly_by_vacc %>% filter(first_presc_yearmonth >= as.Date("2018-01-01")),
       aes(x = first_presc_yearmonth, y = rate * 1000, color = vacc_status)) +
  geom_line(size = 1) +
  annotate("rect",
           xmin = as.Date("2020-01-01"), xmax = as.Date("2021-12-31"),
           ymin = -Inf, ymax = Inf, alpha = 0.2, fill = "gray50") +
  scale_y_continuous(name = "Počet prvopředpisů na 1000 osob bez předpisu", limits = c(0, 3)) +
  labs(
    title = "Prvopředpisovost v čase podle vakcinace",
    subtitle = "Šedá = covidové období (vyřazeno z t-testu)",
    x = "Datum",
    color = "Vakcinace"
  ) +
  theme_minimal()

monthly_rates <- monthly_by_vacc %>%
  filter(period_group %in% c("Pre-COVID", "Post-COVID")) %>%
  mutate(
    vaccinated = if_else(vacc_status == "Vakcinace", 1, 0),
    post_covid = if_else(first_presc_yearmonth >= as.Date("2022-01-01"), 1, 0),
    interaction = vaccinated * post_covid)

did_model <- lm(log(rate * 1000) ~ vaccinated * post_covid, data = monthly_rates)

cat("```\n", paste(capture.output(summary(did_model)), collapse = "\n"), "\n```")
 
Call:
lm(formula = log(rate * 1000) ~ vaccinated * post_covid, data = monthly_rates)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.28041 -0.09699 -0.01711  0.09465  0.48581 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)            0.54788    0.01810  30.266  < 2e-16 ***
vaccinated             0.10123    0.02560   3.954 0.000114 ***
post_covid             0.15113    0.03387   4.463  1.5e-05 ***
vaccinated:post_covid  0.02799    0.04789   0.585 0.559682    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1402 on 164 degrees of freedom
Multiple R-squared:  0.3091,    Adjusted R-squared:  0.2965 
F-statistic: 24.46 on 3 and 164 DF,  p-value: 3.923e-13
 

5.1.4 Multipředpisovost

5.1.4.1 Míra retence (meziroční multipředpisovost)

presc_id_counts <- dbGetQuery(con,
                          "
                          SELECT
                            c.client_id,
                            c.n_vacc,
                            c.n_presc,
                            c.sex,
                            c.birthyear,
                            2021 - c.birthyear AS age_2021,
                            p.date as presc_date,
                            p.presc_order,
                            'CPZP' AS ins_company,
                            YEAR(FIRST(p.date) OVER (PARTITION BY c.client_id ORDER BY p.presc_order ASC)) AS first_presc_year
                          FROM
                            cpzp_clients c
                          LEFT JOIN cpzp_ins_start_end i ON c.client_id = i.client_id
                          LEFT JOIN cpzp_prescriptions p ON c.client_id = p.client_id
                          WHERE
                            i.ins_start_date < '2015-01-01' AND i.ins_end_date = '9999-12-31'
                            AND 2021 - c.birthyear BETWEEN 16 AND 29
                            
                          UNION
                          
                          SELECT
                            c.client_id,
                            c.n_vacc,
                            c.n_presc,
                            c.sex,
                            c.birthyear,
                            2021 - c.birthyear AS age_2021,
                            p.date as presc_date,
                            p.presc_order,
                            'OZP' AS ins_company,
                            YEAR(FIRST(p.date) OVER (PARTITION BY c.client_id ORDER BY p.presc_order ASC)) AS first_presc_year
                          FROM
                            ozp_clients c
                          LEFT JOIN ozp_ins_start_end i ON c.client_id = i.client_id
                          LEFT JOIN ozp_prescriptions p ON c.client_id = p.client_id
                          WHERE
                            i.ins_start_date < '2015-01-01' AND i.ins_end_date IS NULL
                            AND 2021 - c.birthyear BETWEEN 16 AND 29
                          ")%>%
  mutate(presc_year = as.integer(format(as.Date(presc_date), "%Y")),
         vacc_status = ifelse(n_vacc > 0, "Vakcinace", "Bez vakcinace"))

monthly_unique_clients <- presc_id_counts %>%
  mutate(yearmonth = floor_date(presc_date, unit = "month")) %>%
  distinct(client_id, yearmonth, vacc_status) %>%
  count(yearmonth, vacc_status, name = "n_unique_clients")

presc_years <- presc_id_counts %>%
  mutate(
    presc_year = year(presc_date)
  ) %>%
  distinct(client_id, first_presc_year, presc_year, vacc_status) %>%
  filter(first_presc_year >= 2017, first_presc_year <= 2023)


retention_data <- presc_years %>%
  mutate(year_offset = presc_year - first_presc_year) %>%
  filter(year_offset >= 0)


retention_summary <- retention_data %>%
  group_by(first_presc_year, year_offset, vacc_status) %>%
  summarise(n_clients = n_distinct(client_id), .groups = "drop")


cohort_sizes <- retention_data %>%
  filter(year_offset == 0) %>%
  group_by(first_presc_year, vacc_status) %>%
  summarise(cohort_size = n_distinct(client_id), .groups = "drop")


retention_summary <- retention_summary %>%
  left_join(cohort_sizes, by = c("first_presc_year", "vacc_status")) %>%
  mutate(retention_rate = n_clients / cohort_size)

retention_summary_filtered <- retention_summary %>%
  filter(year_offset > 0)

Retence je všeobecně vyšší než u mladší populace a u očkované skupiny je zřejmá vyšší míra udržení “do dalších let”. Reálné zhodnocení vlivu covidu však bude možné až v budoucnosti.

ggplot(retention_summary_filtered, aes(x = year_offset, y = retention_rate, color = as.factor(first_presc_year))) +
  geom_line(size = 1.2) +
  geom_point(size = 2) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_color_viridis_d(name = "Rok prvopředpisu") +
  labs(
    title = "Retence (meziroční předpisovost) podle vakcinace",
    x = "Počet let od prvopředpisu",
    y = "Míra multipředpisovosti"
  ) +
  facet_wrap(~ vacc_status) +
  theme_minimal(base_size = 13)

5.1.4.2 Roční multipředpisovost

multi_presc_yearly <- multi_presc %>%
  mutate(year = year(date)) %>%
  filter(year > 2017)


multi_clients_yearly <- multi_presc_yearly %>%
  group_by(year, vacc_status, client_id) %>%
  summarise(presc_count = n(), .groups = "drop") %>%
  mutate(is_multi = presc_count > 1)


multi_summary_yearly <- multi_clients_yearly %>%
  group_by(year, vacc_status) %>%
  summarise(
    total_clients = n(),
    multi_clients = sum(is_multi),
    multi_rate = (multi_clients / total_clients) * 1000,
    .groups = "drop"
  )


multi_summary_yearly <- multi_summary_yearly %>%
  mutate(
    period_group = case_when(
      year < 2020 ~ "Pre-COVID",
      year > 2021 ~ "Post-COVID",
      TRUE ~ "COVID-Shock"
    )
  )


multi_summary_clean <- multi_summary_yearly %>%
  filter(period_group %in% c("Pre-COVID", "Post-COVID"))


t_test_results <- multi_summary_clean %>%
  group_by(vacc_status) %>%
  summarise(
    p_value = t.test(multi_rate ~ period_group)$p.value,
    mean_before = mean(multi_rate[period_group == "Pre-COVID"]),
    mean_after = mean(multi_rate[period_group == "Post-COVID"]),
    rate = mean_after / mean_before,
    .groups = "drop"
  )

V obou skupinách došlo k nárůstu od původního trendu s rozdílem toho, že v očkované populaci došlo nejprve k poklesu během covidu. Rozdíly mezi skupinami není možné vyhodnotit jako významné.

t_test_results %>%
  transmute(Vakcinace = vacc_status, p_value, `Průměr před` = mean_before, `Průměr po` = mean_after, `Poměr` = rate) %>%
  kable(format = "html", caption = "Porovnání roční multipředpisovosti v čase") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání roční multipředpisovosti v čase
Vakcinace p_value Průměr před Průměr po Poměr
Bez vakcinace 0.0716040 271.5376 311.866 1.148518
Vakcinace 0.1980032 259.4338 288.016 1.110172
ggplot(multi_summary_yearly, aes(x = year, y = multi_rate, color = vacc_status)) +
  geom_line(size = 0.6, alpha = 0.4) +
  stat_smooth(method = "loess", se = FALSE, size = 1.2) +
  annotate("rect",
           xmin = 2020, xmax = 2021,
           ymin = -Inf, ymax = Inf,
           alpha = 0.2, fill = "gray50") +
  scale_y_continuous(name = "Roční multipředpisovost na 1000 předepsaných klientů") +
  scale_x_continuous(breaks = unique(multi_summary_yearly$year)) +
  labs(
    title = "Roční multipředpisovost podle vakcinace",
    subtitle = "Šedá = covidové období (vyřazeno z t-testu)",
    x = "Rok",
    color = "Vakcinace"
  ) +
  theme_minimal()

did_data <- multi_summary_clean %>%
  mutate(
    vaccinated = if_else(vacc_status == "Vakcinace", 1, 0),
    post_covid = if_else(year >= 2022, 1, 0),
    interaction = vaccinated * post_covid
  )

did_model <- lm(log(multi_rate) ~ vaccinated * post_covid, data = did_data)
cat("```\n", paste(capture.output(summary(did_model)), collapse = "\n"), "\n```")
 
Call:
lm(formula = log(multi_rate) ~ vaccinated * post_covid, data = did_data)

Residuals:
       1        2        3        4        5        6        7        8 
-0.02632 -0.04354  0.02632  0.04354 -0.01038 -0.01856  0.01038  0.01856 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)            5.60375    0.02757 203.247 3.52e-09 ***
vaccinated            -0.04620    0.03899  -1.185   0.3016    
post_covid             0.13877    0.03899   3.559   0.0236 *  
vaccinated:post_covid -0.03348    0.05514  -0.607   0.5766    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.03899 on 4 degrees of freedom
Multiple R-squared:  0.8629,    Adjusted R-squared:   0.76 
F-statistic: 8.389 on 3 and 4 DF,  p-value: 0.03361
 

5.1.4.3 Měsíční multipředpisovost

Změny v měsíční předpisovosti v zásadě jen potvrzují předchozí sledování s tím rozdílem, že u očkované skupiny došlo jen k mírnému růstu. Ovšem v tomto případě je rozdíl mezi rozdíly statisticky signifikantní (interakce vaccinated:post_covid).

t_test_results %>%
  transmute(Vakcinace = vacc_status, p_value, `Průměr před` = mean_before, `Průměr po` = mean_after, `Poměr` = rate) %>%
  kable(format = "html", caption = "Porovnání měsíční multipředpisovosti v čase") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání měsíční multipředpisovosti v čase
Vakcinace p_value Průměr před Průměr po Poměr
Bez vakcinace 0.0040403 149.3460 168.9455 1.131235
Vakcinace 0.7959155 147.3262 148.4100 1.007356
ggplot(multi_summary_monthly, aes(x = year_month, y = multi_rate, color = vacc_status)) +
  geom_line(size = 0.5, alpha = 0.4) + 
  stat_smooth(method = "loess", se = FALSE, size = 1.2) +  
  annotate("rect",
           xmin = as.Date("2020-01-01"), xmax = as.Date("2021-12-31"),
           ymin = -Inf, ymax = Inf, alpha = 0.2, fill = "gray50") +
  scale_y_continuous(name = "Měsíční multipředpisovost na 1000 předepsaných klientů") +
  labs(
    title = "Měsíční multipředpisovost podle vakcinačního statusu",
    subtitle = "Šedá = covidové období (vyřazeno z t-testu)",
    x = "Datum",
    color = "Vakcinace"
  ) +
  theme_minimal()

did_data <- multi_summary_clean %>%
  mutate(
    vaccinated = if_else(vacc_status == "Vakcinace", 1, 0),
    post_covid = if_else(year_month >= as.Date("2022-01-01"), 1, 0),
    interaction = vaccinated * post_covid
  )

did_model <- lm(log(multi_rate) ~ vaccinated * post_covid, data = did_data)
cat("```\n", paste(capture.output(summary(did_model)), collapse = "\n"), "\n```")
 
Call:
lm(formula = log(multi_rate) ~ vaccinated * post_covid, data = did_data)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.266633 -0.081329  0.002076  0.091950  0.284863 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)            4.994406   0.024953 200.152  < 2e-16 ***
vaccinated            -0.007215   0.035289  -0.204 0.838443    
post_covid             0.127420   0.035289   3.611 0.000497 ***
vaccinated:post_covid -0.118159   0.049906  -2.368 0.019997 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1222 on 92 degrees of freedom
Multiple R-squared:  0.1798,    Adjusted R-squared:  0.153 
F-statistic: 6.722 on 3 and 92 DF,  p-value: 0.0003761
 

5.1.5 Způsob podání

Vývoj zde jen potvrzuje trend nastolený v předchozí sledované skupině a tedy, že s rostoucím věkem se preferují injekce / infuze na úkor tablet.

prescriptions <- prescriptions %>%
  mutate(
    drug_form_std = case_when(
      is.na(drug_form) ~ NA_character_,
      str_detect(drug_form, "tableta|Tableta") ~ "Tableta",
      str_detect(drug_form, "tobolka|Tobolka") ~ "Tobolka",
      str_detect(drug_form, "injekční|Injekční|infuzní|Infuzní") ~ "Injekce / Infuze",
      str_detect(drug_form, "č[íi]pek|Čípek") ~ "Čípek",
      str_detect(drug_form, "perorální roztok|Perorální roztok") ~ "Perorální roztok",
      TRUE ~ "Other"
    )
  )

prescriptions <- prescriptions %>%
  mutate(
    year = year(presc_date)
  )

ggplot(prescriptions %>% filter(!is.na(presc_date), !is.na(drug_form)), aes(x = factor(year), fill = drug_form_std)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(
    title = "Poměr způsobu podání v čase",
    x = "Rok",
    y = "%",
    fill = "Způsob podání"
  ) +
  theme_minimal()

client_yearly_counts <- prescriptions %>%
  filter(!is.na(presc_date), !is.na(drug_form)) %>%
  group_by(client_id, year, drug_form_std) %>%
  summarise(presc_count = n(), .groups = "drop")

average_freq_per_year <- client_yearly_counts %>%
  group_by(year, drug_form_std) %>%
  summarise(avg_presc_per_client = mean(presc_count), .groups = "drop")

ggplot(average_freq_per_year, aes(x = factor(year), y = avg_presc_per_client, fill = drug_form_std)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(
    title = "Počet předpisů za rok na klienta dle způsobu podání",
    x = "Rok",
    y = "Průměrný počet předpisů na klienta",
    fill = "Způsob podání"
  ) +
  theme_minimal()

Při sledování poměru předpisů je zajímavé sledovat výkyv neočkované populace ve druhé polovině roku 2021. K tomuto výkyvu dochází u všech dospělých neočkovaných skupin a je způsoben jak snížením celkového množství injekcí, tak zvýšením celkového množství tablet. Může se jednat o vlivy covidových nařízení a částečné nepřístupnosti ambulantní léčby pro neočkovanou populaci.

filtered_data <- prescriptions %>%
  mutate(
    drug_form_std = ifelse(is.na(drug_form_std), "Unknown", drug_form_std)
  ) %>%
  filter(
    drug_form_std %in% c("Injekce / Infuze", "Tableta"),
    !is.na(vaccinated),
    !is.na(presc_yearmonth)
  )

form_vax_counts <- filtered_data %>%
  group_by(presc_yearmonth, vaccinated, drug_form_std) %>%
  summarise(n = n(), .groups = "drop")

month_totals <- filtered_data %>%
  group_by(presc_yearmonth, vaccinated) %>%
  summarise(total = n(), .groups = "drop")

form_vax_ratios <- form_vax_counts %>%
  left_join(month_totals, by = c("presc_yearmonth", "vaccinated")) %>%
  mutate(ratio = n / total)

form_vax_ratios <- form_vax_ratios %>%
  mutate(group = paste(drug_form_std, ifelse(vaccinated == "Vakcinace", "Vakcinace", "Bez vakcinace"), sep = " - "))

ggplot(form_vax_ratios, aes(x = presc_yearmonth, y = ratio, color = group)) +
  geom_line(size = 0.6) +                   
  geom_smooth(se = FALSE, span = 0.3) +   
  scale_y_continuous(labels = percent_format(accuracy = 1)) +
  scale_x_date(date_labels = "%Y",         
               date_breaks = "1 year",
               expand = c(0.01, 0)) +
  labs(
    title = "Poměr počtu předpisů podle způsobu podání a vakcinace",
    x = "Datum",
    y = "Poměr počtu předpisů",
    color = "Forma podání + Vakcinace"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0, hjust = 0.5))

5.2 Shrnutí

V rámci dlouhodobého vývoje jsou opět všechna rizika větší v době po covidu oproti době před covidem. Průběh očkovaných a neočkovaných je však srovnatelný. Jedinou výjimkou je měsíční multipředpisovost, kde očkovaní podléhají mírně nižšímu riziku. Poměry způsobu podání v zásadě kopírují trend nastavený v mladší kategorii - a tedy nárůst injekcí / infuzí na úkor tablet s jedinou výjimkou ve druhé polovině roku 2021, kde výrazně rostou tablety a klesají injekce / infuze.

Krátkodobá analýza je sloučená se skupinou 30-39 z důvodu nízké relevance výsledků.

6 Věková kategorie 30-39

6.1 Dlouhodobý efekt

6.1.1 Distribuce věku a pohlaví

Demografické rozložení je v tomto případě mírně vyváženější než u mladší populace.

unique_clients <- prescriptions %>%
  filter(!is.na(age_2021), !is.na(sex), !is.na(vaccinated)) %>%
  distinct(client_id, .keep_all = TRUE)  

age_data <- unique_clients %>%
  group_by(vaccinated, age_2021) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(vaccinated) %>%
  mutate(prop = count / sum(count),
         variable = "Age",
         value = as.character(age_2021)) %>%
  ungroup()

sex_data <- unique_clients %>%
  group_by(vaccinated, sex) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(vaccinated) %>%
  mutate(prop = count / sum(count),
         variable = "Sex",
         value = sex) %>%
  ungroup()
ggplot(unique_clients, aes(x = age_2021, fill = factor(vaccinated))) +
  geom_histogram(aes(y = after_stat(density)), position = "identity", alpha = 0.4, binwidth = 1) +
  scale_fill_brewer(palette = "Set1", name = "Vaccinated") +
  scale_x_continuous(breaks = seq(floor(min(unique_clients$age)), ceiling(max(unique_clients$age)), by = 1)) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    title = "Věkové rozložení podle statusu vakcinace",
    x = "Věk", y = "%",
    fill = "Vakcinace"
  ) +
  theme_minimal()

ggplot(sex_data, aes(x = sex, y = prop, color = vaccinated, group = vaccinated)) +
  geom_line(size = 1) +
  geom_point(size = 2) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    title = "Pohlaví podle vakcinace",
    x = "Pohlaví",
    y = "%",
    color = "Vakcinace"
  ) +
  theme_minimal()

6.1.2 Předpověď vs realita

prescriptions_clean <- prescriptions %>%
  filter(!is.na(presc_date)) %>%
  mutate(presc_date = as.Date(presc_date)) %>%
  arrange(presc_date) %>%
  filter(presc_date >= as.Date("2015-01-01") & presc_date <= as.Date("2023-12-31")) %>%
  mutate(count = 1)  


results_vacc <- prescriptions_clean %>%
  filter(!is.na(vaccinated)) %>%
  group_by(vaccinated) %>%
  group_split() %>%
  map(~ {
    group_val <- unique(.x$vaccinated)
    analyze_group(df = filter(prescriptions_clean, vaccinated == group_val), 
                  group_id = group_val, 
                  group_type = "Vaccination")
  }) %>%
  compact() %>%
  bind_rows()

results_sex <- prescriptions_clean %>%
  filter(!is.na(sex)) %>%
  group_by(sex) %>%
  group_split() %>%
  map(~ {
    group_val <- unique(.x$sex)
    analyze_group(df = filter(prescriptions_clean, sex == group_val), 
                  group_id = group_val, 
                  group_type = "Sex")
  }) %>%
  compact() %>%
  bind_rows()

results_all <- bind_rows(results_vacc, results_sex) %>%
  arrange(p_value)

plot_data_all <- results_all %>%
  unnest(plot_data, names_sep = "_")

Rozdíly očekávaného úhrnu předpisů jsou napříč všemi kategoriemi téměř identické.

results_all %>%
  mutate(Skupina = group, RMSE = rmse, rRMSE = rrmse) %>%
  select(
    Skupina, RMSE, rRMSE, t_stat, p_value
  ) %>%
  kable(format = "html", caption = "Porovnání předpisovosti v čase (ETS)") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání předpisovosti v čase (ETS)
Skupina RMSE rRMSE t_stat p_value
M 78.20163 0.1367719 -3.374315 0.0014910
Bez vakcinace 59.18631 0.1257455 -2.990672 0.0044209
F 118.90336 0.1315060 -2.603776 0.0123016
Vakcinace 127.09959 0.1264358 -2.449810 0.0180735
ggplot(plot_data_all, aes(x = plot_data_month)) +
  geom_line(aes(y = plot_data_prescriptions), color = "blue") +
  geom_line(aes(y = plot_data_predicted), color = "red", linetype = "dashed") +
  geom_vline(xintercept = as.Date("2020-01-01"), linetype = "dotted") +
  facet_wrap(~ paste(plot_data_group_type, plot_data_group, sep = ": "), scales = "free_y") +
  labs(
    title = "Realita vs Předpověď (ETS)",
    subtitle = "Modrá = Realita, Červená = Předpověď",
    x = "Datum", y = "Počet předpisů"
  ) +
  theme_minimal()

6.1.3 Prvopředpisovost

first_presc <- dbGetQuery(con,
                          "
                          SELECT
                            c.client_id,
                            c.n_vacc,
                            c.n_presc,
                            c.sex,
                            c.birthyear,
                            FIRST(p.date) OVER (PARTITION BY c.client_id ORDER BY p.date ASC) AS first_presc_date,
                            'CPZP' AS ins_company
                          FROM
                            cpzp_clients c
                          LEFT JOIN cpzp_ins_start_end i ON c.client_id = i.client_id
                          LEFT JOIN cpzp_prescriptions p ON c.client_id = p.client_id AND p.presc_order = 1
                          WHERE
                            i.ins_start_date < '2015-01-01' AND i.ins_end_date = '9999-12-31'
                            AND 2021 - c.birthyear BETWEEN 30 AND 39
                            
                          UNION
                          
                          SELECT
                            c.client_id,
                            c.n_vacc,
                            c.n_presc,
                            c.sex,
                            c.birthyear,
                            FIRST(p.date) OVER (PARTITION BY c.client_id ORDER BY p.date ASC) AS first_presc_date,
                            'OZP' AS ins_company
                          FROM
                            ozp_clients c
                          LEFT JOIN ozp_ins_start_end i ON c.client_id = i.client_id
                          LEFT JOIN ozp_prescriptions p ON c.client_id = p.client_id AND p.presc_order = 1
                          WHERE
                            i.ins_start_date < '2015-01-01' AND i.ins_end_date IS NULL
                            AND 2021 - c.birthyear BETWEEN 30 AND 39
                          ") %>%
  mutate(first_presc_yearmonth = floor_date(first_presc_date, unit = "month"),
         sex = ifelse(sex == "Z", "F", sex),
         vacc_status = ifelse(n_vacc > 0, "Vakcinace", "Bez vakcinace"),
         age_2021 = 2021 - birthyear)

monthly_first_presc <- first_presc %>%
  count(first_presc_yearmonth, name = "new_prescriptions") %>%
  arrange(first_presc_yearmonth) %>%
  mutate(
    cum_new_presc = cumsum(new_prescriptions),  
    total_population = nrow(first_presc),
    at_risk_population = total_population - lag(cum_new_presc, default = 0),
    rate = new_prescriptions / at_risk_population
  ) %>%
  mutate(
    period_group = case_when(
      first_presc_yearmonth < as.Date("2020-01-01") ~ "Pre-COVID",
      first_presc_yearmonth > as.Date("2021-12-31") ~ "Post-COVID",
      TRUE ~ "COVID-Shock"
    )
  )

monthly_clean <- monthly_first_presc %>%
  filter(first_presc_yearmonth >= as.Date("2018-01-01")) %>%
  filter(period_group %in% c("Pre-COVID", "Post-COVID"))

ttest <- t.test(rate ~ period_group, data = monthly_clean)

Stejně jako u jiných věkových skupin zde došlo k významnému nárůstu prvopředpisovosti.

ggplot(monthly_first_presc %>% filter(first_presc_yearmonth >= as.Date("2018-01-01")),
       aes(x = first_presc_yearmonth, y = rate * 1000)) +
  geom_line(color = "steelblue", size = 1) +
  annotate("rect",
           xmin = as.Date("2020-01-01"), xmax = as.Date("2021-12-31"),
           ymin = -Inf, ymax = Inf, alpha = 0.2, fill = "gray50") +
  scale_y_continuous(name = "Počet prvopředpisů na 1000 osob bez předpisu", limits = c(0, 5)) +
  labs(
    title = "Prvopředpisovost v čase",
    subtitle = "Šedá = covidové období (vyřazeno z t-testu)",
    x = "Datum"
  ) +
  theme_minimal()

cat("```\n", paste(capture.output(ttest), collapse = "\n"), "\n```")
 
    Welch Two Sample t-test

data:  rate by period_group
t = 4.3496, df = 45.993, p-value = 7.499e-05
alternative hypothesis: true difference in means between group Post-COVID and group Pre-COVID is not equal to 0
95 percent confidence interval:
 0.0001565338 0.0004262142
sample estimates:
mean in group Post-COVID  mean in group Pre-COVID 
             0.003192046              0.002900672 
 
total_by_vacc <- first_presc %>%
  count(vacc_status, name = "total_population")

monthly_by_vacc <- first_presc %>%
  count(vacc_status, first_presc_yearmonth, name = "new_prescriptions") %>%
  arrange(vacc_status, first_presc_yearmonth) %>%
  group_by(vacc_status) %>%
  mutate(
    cum_prescriptions = cumsum(new_prescriptions),
    total_population = total_by_vacc$total_population[match(vacc_status, total_by_vacc$vacc_status)],
    at_risk_population = total_population - lag(cum_prescriptions, default = 0),
    rate = new_prescriptions / at_risk_population
  ) %>%
  ungroup() %>%
  mutate(
    period_group = case_when(
      first_presc_yearmonth < as.Date("2020-01-01") ~ "Pre-COVID",
      first_presc_yearmonth > as.Date("2021-12-31") ~ "Post-COVID",
      TRUE ~ "COVID-Shock"
    )
  )
  
monthly_by_vacc_clean <- monthly_by_vacc %>%
  filter(first_presc_yearmonth >= as.Date("2018-01-01")) %>%
  filter(period_group %in% c("Pre-COVID", "Post-COVID"))

first_presc_table <- monthly_by_vacc_clean %>%
  group_by(vacc_status) %>%
  summarise(
    p_value = t.test(rate ~ period_group)$p.value,
    mean_before = mean(rate[period_group == "Pre-COVID"]),
    mean_after = mean(rate[period_group == "Post-COVID"]),
    rate = mean_after / mean_before
  )

K tomuto nárůstu došlo napříč vakcinovanou i nevakcinovanou populací a rozdíly mezi populacemi nejsou významné.

first_presc_table %>%
  transmute(Vakcinace = vacc_status, p_value, `Průměr před` = mean_before * 1000, `Průměr po` = mean_after * 1000, `Poměr` = rate) %>%
  kable(format = "html", caption = "Porovnání prvopředpisovosti v čase") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání prvopředpisovosti v čase
Vakcinace p_value Průměr před Průměr po Poměr
Bez vakcinace 0.0197657 2.625809 2.816226 1.072518
Vakcinace 0.0000141 3.055007 3.407820 1.115487
ggplot(monthly_by_vacc %>% filter(first_presc_yearmonth >= as.Date("2018-01-01")),
       aes(x = first_presc_yearmonth, y = rate * 1000, color = vacc_status)) +
  geom_line(size = 1) +
  annotate("rect",
           xmin = as.Date("2020-01-01"), xmax = as.Date("2021-12-31"),
           ymin = -Inf, ymax = Inf, alpha = 0.2, fill = "gray50") +
  scale_y_continuous(name = "Počet prvopředpisů na 1000 osob bez předpisu", limits = c(0, 5)) +
  labs(
    title = "Prvopředpisovost v čase podle vakcinace",
    subtitle = "Šedá = covidové období (vyřazeno z t-testu)",
    x = "Datum",
    color = "Vakcinace"
  ) +
  theme_minimal()

monthly_rates <- monthly_by_vacc %>%
  filter(period_group %in% c("Pre-COVID", "Post-COVID")) %>%
  mutate(
    vaccinated = if_else(vacc_status == "Vakcinace", 1, 0),
    post_covid = if_else(first_presc_yearmonth >= as.Date("2022-01-01"), 1, 0),
    interaction = vaccinated * post_covid)

did_model <- lm(log(rate * 1000) ~ vaccinated * post_covid, data = monthly_rates)

cat("```\n", paste(capture.output(summary(did_model)), collapse = "\n"), "\n```")
 
Call:
lm(formula = log(rate * 1000) ~ vaccinated * post_covid, data = monthly_rates)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.37499 -0.07649 -0.01110  0.06480  0.48692 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)           1.026739   0.017329  59.249  < 2e-16 ***
vaccinated            0.133896   0.024507   5.464  1.7e-07 ***
post_covid            0.003442   0.032420   0.106    0.916    
vaccinated:post_covid 0.059457   0.045849   1.297    0.197    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1342 on 164 degrees of freedom
Multiple R-squared:  0.2574,    Adjusted R-squared:  0.2438 
F-statistic: 18.95 on 3 and 164 DF,  p-value: 1.338e-10
 

6.1.4 Multipředpisovost

6.1.4.1 Míra retence (meziroční multipředpisovost)

Zde je opět potvrzený trend vyšší míry oproti mladší populaci a také vyšší retence očkovaných jedinců.

presc_id_counts <- dbGetQuery(con,
                          "
                          SELECT
                            c.client_id,
                            c.n_vacc,
                            c.n_presc,
                            c.sex,
                            c.birthyear,
                            2021 - c.birthyear AS age_2021,
                            p.date as presc_date,
                            p.presc_order,
                            'CPZP' AS ins_company,
                            YEAR(FIRST(p.date) OVER (PARTITION BY c.client_id ORDER BY p.presc_order ASC)) AS first_presc_year
                          FROM
                            cpzp_clients c
                          LEFT JOIN cpzp_ins_start_end i ON c.client_id = i.client_id
                          LEFT JOIN cpzp_prescriptions p ON c.client_id = p.client_id
                          WHERE
                            i.ins_start_date < '2015-01-01' AND i.ins_end_date = '9999-12-31'
                            AND 2021 - c.birthyear BETWEEN 30 AND 39
                            
                          UNION
                          
                          SELECT
                            c.client_id,
                            c.n_vacc,
                            c.n_presc,
                            c.sex,
                            c.birthyear,
                            2021 - c.birthyear AS age_2021,
                            p.date as presc_date,
                            p.presc_order,
                            'OZP' AS ins_company,
                            YEAR(FIRST(p.date) OVER (PARTITION BY c.client_id ORDER BY p.presc_order ASC)) AS first_presc_year
                          FROM
                            ozp_clients c
                          LEFT JOIN ozp_ins_start_end i ON c.client_id = i.client_id
                          LEFT JOIN ozp_prescriptions p ON c.client_id = p.client_id
                          WHERE
                            i.ins_start_date < '2015-01-01' AND i.ins_end_date IS NULL
                            AND 2021 - c.birthyear BETWEEN 30 AND 39
                          ")%>%
  mutate(presc_year = as.integer(format(as.Date(presc_date), "%Y")),
         vacc_status = ifelse(n_vacc > 0, "Vakcinace", "Bez vakcinace"))

monthly_unique_clients <- presc_id_counts %>%
  mutate(yearmonth = floor_date(presc_date, unit = "month")) %>%
  distinct(client_id, yearmonth, vacc_status) %>%
  count(yearmonth, vacc_status, name = "n_unique_clients")

presc_years <- presc_id_counts %>%
  mutate(
    presc_year = year(presc_date)
  ) %>%
  distinct(client_id, first_presc_year, presc_year, vacc_status) %>%
  filter(first_presc_year >= 2017, first_presc_year <= 2023)


retention_data <- presc_years %>%
  mutate(year_offset = presc_year - first_presc_year) %>%
  filter(year_offset >= 0)


retention_summary <- retention_data %>%
  group_by(first_presc_year, year_offset, vacc_status) %>%
  summarise(n_clients = n_distinct(client_id), .groups = "drop")


cohort_sizes <- retention_data %>%
  filter(year_offset == 0) %>%
  group_by(first_presc_year, vacc_status) %>%
  summarise(cohort_size = n_distinct(client_id), .groups = "drop")


retention_summary <- retention_summary %>%
  left_join(cohort_sizes, by = c("first_presc_year", "vacc_status")) %>%
  mutate(retention_rate = n_clients / cohort_size)

retention_summary_filtered <- retention_summary %>%
  filter(year_offset > 0)
ggplot(retention_summary_filtered, aes(x = year_offset, y = retention_rate, color = as.factor(first_presc_year))) +
  geom_line(size = 1.2) +
  geom_point(size = 2) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_color_viridis_d(name = "Rok prvopředpisu") +
  labs(
    title = "Retence (meziroční předpisovost) podle vakcinace",
    x = "Počet let od prvopředpisu",
    y = "Míra multipředpisovosti"
  ) +
  facet_wrap(~ vacc_status) +
  theme_minimal(base_size = 13)

6.1.4.2 Roční multipředpisovost

Během porovnání obou období dochází k mírnému nárůstu. Zde je ale nutné zmínit, že samotný průměr míry multipředpisovosti není vhodnou metrikou a to z důvodu stoupavého trendu až do roku 2020 a poté spíše konstantního průběhu pro očkovanou populaci. Rozdíly mezi skupinami jsou opět nevýznamné.

multi_presc_yearly <- multi_presc %>%
  mutate(year = year(date)) %>%
  filter(year > 2017)


multi_clients_yearly <- multi_presc_yearly %>%
  group_by(year, vacc_status, client_id) %>%
  summarise(presc_count = n(), .groups = "drop") %>%
  mutate(is_multi = presc_count > 1)


multi_summary_yearly <- multi_clients_yearly %>%
  group_by(year, vacc_status) %>%
  summarise(
    total_clients = n(),
    multi_clients = sum(is_multi),
    multi_rate = (multi_clients / total_clients) * 1000,
    .groups = "drop"
  )


multi_summary_yearly <- multi_summary_yearly %>%
  mutate(
    period_group = case_when(
      year < 2020 ~ "Pre-COVID",
      year > 2021 ~ "Post-COVID",
      TRUE ~ "COVID-Shock"
    )
  )


multi_summary_clean <- multi_summary_yearly %>%
  filter(period_group %in% c("Pre-COVID", "Post-COVID"))


t_test_results <- multi_summary_clean %>%
  group_by(vacc_status) %>%
  summarise(
    p_value = t.test(multi_rate ~ period_group)$p.value,
    mean_before = mean(multi_rate[period_group == "Pre-COVID"]),
    mean_after = mean(multi_rate[period_group == "Post-COVID"]),
    rate = mean_after / mean_before,
    .groups = "drop"
  )
t_test_results %>%
  transmute(Vakcinace = vacc_status, p_value, `Průměr před` = mean_before, `Průměr po` = mean_after, `Poměr` = rate) %>%
  kable(format = "html", caption = "Porovnání roční multipředpisovosti v čase") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání roční multipředpisovosti v čase
Vakcinace p_value Průměr před Průměr po Poměr
Bez vakcinace 0.2185289 346.9531 359.5639 1.036347
Vakcinace 0.1433689 336.7882 358.9081 1.065679
ggplot(multi_summary_yearly, aes(x = year, y = multi_rate, color = vacc_status)) +
  geom_line(size = 0.6, alpha = 0.4) +
  stat_smooth(method = "loess", se = FALSE, size = 1.2) +
  annotate("rect",
           xmin = 2020, xmax = 2021,
           ymin = -Inf, ymax = Inf,
           alpha = 0.2, fill = "gray50") +
  scale_y_continuous(name = "Roční multipředpisovost na 1000 předepsaných klientů") +
  scale_x_continuous(breaks = unique(multi_summary_yearly$year)) +
  labs(
    title = "Roční multipředpisovost podle vakcinace",
    subtitle = "Šedá = covidové období (vyřazeno z t-testu)",
    x = "Rok",
    color = "Vakcinace"
  ) +
  theme_minimal()

did_data <- multi_summary_clean %>%
  mutate(
    vaccinated = if_else(vacc_status == "Vakcinace", 1, 0),
    post_covid = if_else(year >= 2022, 1, 0),
    interaction = vaccinated * post_covid
  )

did_model <- lm(log(multi_rate) ~ vaccinated * post_covid, data = did_data)
cat("```\n", paste(capture.output(summary(did_model)), collapse = "\n"), "\n```")
 
Call:
lm(formula = log(multi_rate) ~ vaccinated * post_covid, data = did_data)

Residuals:
        1         2         3         4         5         6         7         8 
-0.002338 -0.018123  0.002338  0.018123  0.013060 -0.005682 -0.013060  0.005682 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)            5.84919    0.01158 504.933 9.23e-11 ***
vaccinated            -0.02990    0.01638  -1.825   0.1421    
post_covid             0.03562    0.01638   2.174   0.0954 .  
vaccinated:post_covid  0.02814    0.02317   1.215   0.2913    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.01638 on 4 degrees of freedom
Multiple R-squared:  0.8446,    Adjusted R-squared:  0.7281 
F-statistic: 7.247 on 3 and 4 DF,  p-value: 0.04286
 

6.1.4.3 Měsíční multipředpisovost

Měsíční multipředpisovost vykazuje mírně klesavou tendenci v období po covidu u nevakcinované populace a stoupaní u vakcinované. Rozdíly uvnitř skupin i mezi skupinami jsou však nevýznamné.

t_test_results %>%
  transmute(Vakcinace = vacc_status, p_value, `Průměr před` = mean_before, `Průměr po` = mean_after, `Poměr` = rate) %>%
  kable(format = "html", caption = "Porovnání měsíční multipředpisovosti v čase") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání měsíční multipředpisovosti v čase
Vakcinace p_value Průměr před Průměr po Poměr
Bez vakcinace 0.6580879 152.3939 154.5434 1.014105
Vakcinace 0.7603065 153.0694 154.0509 1.006412
ggplot(multi_summary_monthly, aes(x = year_month, y = multi_rate, color = vacc_status)) +
  geom_line(size = 0.5, alpha = 0.4) + 
  stat_smooth(method = "loess", se = FALSE, size = 1.2) +  
  annotate("rect",
           xmin = as.Date("2020-01-01"), xmax = as.Date("2021-12-31"),
           ymin = -Inf, ymax = Inf, alpha = 0.2, fill = "gray50") +
  scale_y_continuous(name = "Měsíční multipředpisovost na 1000 předepsaných klientů") +
  labs(
    title = "Měsíční multipředpisovost podle vakcinačního statusu",
    subtitle = "Šedá = covidové období (vyřazeno z t-testu)",
    x = "Datum",
    color = "Vakcinace"
  ) +
  theme_minimal()

did_data <- multi_summary_clean %>%
  mutate(
    vaccinated = if_else(vacc_status == "Vakcinace", 1, 0),
    post_covid = if_else(year_month >= as.Date("2022-01-01"), 1, 0),
    interaction = vaccinated * post_covid
  )

did_model <- lm(log(multi_rate) ~ vaccinated * post_covid, data = did_data)
cat("```\n", paste(capture.output(summary(did_model)), collapse = "\n"), "\n```")
 
Call:
lm(formula = log(multi_rate) ~ vaccinated * post_covid, data = did_data)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.264666 -0.053956  0.005046  0.062361  0.290939 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)            5.019074   0.018752 267.662   <2e-16 ***
vaccinated             0.008983   0.026519   0.339    0.736    
post_covid             0.017669   0.026519   0.666    0.507    
vaccinated:post_covid -0.010698   0.037503  -0.285    0.776    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.09186 on 92 degrees of freedom
Multiple R-squared:  0.005949,  Adjusted R-squared:  -0.02647 
F-statistic: 0.1835 on 3 and 92 DF,  p-value: 0.9073
 

6.1.5 Způsob podání

Poměry způsobu podání se oproti mladším skupinám srovnaly a vykazují zhruba stejné hodnoty.

prescriptions <- prescriptions %>%
  mutate(
    drug_form_std = case_when(
      is.na(drug_form) ~ NA_character_,
      str_detect(drug_form, "tableta|Tableta") ~ "Tableta",
      str_detect(drug_form, "tobolka|Tobolka") ~ "Tobolka",
      str_detect(drug_form, "injekční|Injekční|infuzní|Infuzní") ~ "Injekce / Infuze",
      str_detect(drug_form, "č[íi]pek|Čípek") ~ "Čípek",
      str_detect(drug_form, "perorální roztok|Perorální roztok") ~ "Perorální roztok",
      TRUE ~ "Other"
    )
  )

prescriptions <- prescriptions %>%
  mutate(
    year = year(presc_date)
  )

ggplot(prescriptions %>% filter(!is.na(presc_date), !is.na(drug_form)), aes(x = factor(year), fill = drug_form_std)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(
    title = "Poměr způsobu podání v čase",
    x = "Rok",
    y = "%",
    fill = "Způsob podání"
  ) +
  theme_minimal()

client_yearly_counts <- prescriptions %>%
  filter(!is.na(presc_date), !is.na(drug_form)) %>%
  group_by(client_id, year, drug_form_std) %>%
  summarise(presc_count = n(), .groups = "drop")

average_freq_per_year <- client_yearly_counts %>%
  group_by(year, drug_form_std) %>%
  summarise(avg_presc_per_client = mean(presc_count), .groups = "drop")

ggplot(average_freq_per_year, aes(x = factor(year), y = avg_presc_per_client, fill = drug_form_std)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(
    title = "Počet předpisů za rok na klienta dle způsobu podání",
    x = "Rok",
    y = "Průměrný počet předpisů na klienta",
    fill = "Způsob podání"
  ) +
  theme_minimal()

Při vykreslení poměru injekcí a infuzí na tablety je zde opět zajímavý nárůst ve prospěch tablet ve druhé polovině roku 2021.

filtered_data <- prescriptions %>%
  mutate(
    drug_form_std = ifelse(is.na(drug_form_std), "Unknown", drug_form_std)
  ) %>%
  filter(
    drug_form_std %in% c("Injekce / Infuze", "Tableta"),
    !is.na(vaccinated),
    !is.na(presc_yearmonth)
  )

form_vax_counts <- filtered_data %>%
  group_by(presc_yearmonth, vaccinated, drug_form_std) %>%
  summarise(n = n(), .groups = "drop")

month_totals <- filtered_data %>%
  group_by(presc_yearmonth, vaccinated) %>%
  summarise(total = n(), .groups = "drop")

form_vax_ratios <- form_vax_counts %>%
  left_join(month_totals, by = c("presc_yearmonth", "vaccinated")) %>%
  mutate(ratio = n / total)

form_vax_ratios <- form_vax_ratios %>%
  mutate(group = paste(drug_form_std, ifelse(vaccinated == "Vakcinace", "Vakcinace", "Bez vakcinace"), sep = " - "))

ggplot(form_vax_ratios, aes(x = presc_yearmonth, y = ratio, color = group)) +
  geom_line(size = 0.6) +                   
  geom_smooth(se = FALSE, span = 0.3) +   
  scale_y_continuous(labels = percent_format(accuracy = 1)) +
  scale_x_date(date_labels = "%Y",         
               date_breaks = "1 year",
               expand = c(0.01, 0)) +
  labs(
    title = "Poměr počtu předpisů podle způsobu podání a vakcinace",
    x = "Datum",
    y = "Poměr počtu předpisů",
    color = "Forma podání + Vakcinace"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0, hjust = 0.5))

6.2 Shrnutí

Zatímco prvopředpisovost významně stoupla, tak multipředpisovost zůstává v zásadě na stejné úrovni a v některých případech dokonce i klesá. Rozdíly mezi očkovanými a neočkovanými jedinci jsou zanedbatelné.

Trendy způsobu podání jsou v této skupině oproti mladším populacím ustálené s jediným výkyvem ve druhé polovině roku 2022.

7 Krátkodobý efekt 16-39

7.1 Prvopředpisy

7.1.1 SCCS vakcinovaných

Všeobecně je riziko prvopředpisu po vakcinaci vyšší oproti baseline před očkováním. Z hlediska statistické významnosti lze ale s určitou dávkou jistoty takto zhodnotit pouze období 0-29, 90-179 a 270-365.

combined_vaccinated_prescribed <- dbGetQuery(con,
                                    "
                                    SELECT
                                    v.client_id,
                                    v.date AS vaccination_date,
                                    p.date AS event_date,
                                    'CPZP' as ins_company
                                    FROM cpzp_vaccinations v
                                    LEFT JOIN cpzp_clients c ON v.client_id = c.client_id
                                    LEFT JOIN cpzp_ins_start_end s ON c.client_id = s.client_id
                                    INNER JOIN cpzp_prescriptions p ON c.client_id = p.client_id
                                                                    AND p.presc_order = 1
                                                                    AND p.date BETWEEN v.date - INTERVAL '1 year' AND v.date + INTERVAL '1 year'
                                    WHERE v.event_order = 1
                                    AND s.ins_start_date < v.date - INTERVAL '1 year'
                                    AND s.ins_end_date > v.date + INTERVAL '1 year'
                                    AND v.date BETWEEN DATE '2021-06-15' - INTERVAL '7 days'
                                                       AND DATE '2021-06-15' + INTERVAL '7 days'
                                    AND 2021 - c.birthyear BETWEEN 16 AND 39
                                    
                                    UNION
                                    
                                    SELECT
                                    v.client_id,
                                    v.date AS vaccination_date,
                                    p.date AS event_date,
                                    'OZP' AS ins_company
                                    FROM OZP_vaccinations v
                                    LEFT JOIN ozp_clients c ON v.client_id = c.client_id
                                    LEFT JOIN ozp_ins_start_end s ON c.client_id = s.client_id
                                    INNER JOIN ozp_prescriptions p ON c.client_id = p.client_id
                                                                    AND p.presc_order = 1
                                                                    AND p.date BETWEEN v.date - INTERVAL '1 year' AND v.date + INTERVAL '1 year'
                                    WHERE v.event_order = 1
                                    AND s.ins_start_date < v.date - INTERVAL '1 year'
                                    AND s.ins_end_date > v.date + INTERVAL '1 year'
                                    AND v.date BETWEEN DATE '2021-06-15' - INTERVAL '7 days'
                                                       AND DATE '2021-06-15' + INTERVAL '7 days'
                                    AND 2021 - c.birthyear BETWEEN 16 AND 39
                                    ") %>%
  mutate(client_id = row_number()) %>%
  transmute(
    person_id = client_id,
    exposure_start_date = vaccination_date,
    event_date,
    observation_start_date = vaccination_date - days(365),
    observation_end_date = vaccination_date + days(365)
  )

sccs_vaccinated <- combined_vaccinated_prescribed %>%
  mutate(
    event_day = as.integer(event_date - observation_start_date),
    exposure_start_day = as.integer(exposure_start_date - observation_start_date),
    observation_start_day = 0,
    observation_end_day = as.integer(observation_end_date - observation_start_date)
  )

sccs_result_vaccinated <- standardsccs(
  formula = event_day ~ exposure_start_day,
  indiv = person_id,
  astart = observation_start_day,
  aend = observation_end_day,
  aevent = event_day,
  adrug = exposure_start_day,
  aedrug = exposure_start_day + 365,
  expogrp = list(c(0, 30, 90, 180, 270)),
  data = sccs_vaccinated
)

cat("```\n", paste(capture.output(print(sccs_result_vaccinated)), collapse = "\n"), "\n```")
 Call:
coxph(formula = Surv(rep(1, 12816L), event) ~ exposure_start_day + 
    strata(indivL) + offset(log(interval)), data = chopdat, method = "exact")

  n= 12816, number of events= 2136 

                       coef exp(coef) se(coef)     z Pr(>|z|)  
exposure_start_day1 0.25395   1.28910  0.10263 2.474   0.0133 *
exposure_start_day2 0.11084   1.11722  0.08064 1.374   0.1693  
exposure_start_day3 0.14328   1.15405  0.06749 2.123   0.0338 *
exposure_start_day4 0.10718   1.11313  0.06845 1.566   0.1174  
exposure_start_day5 0.15386   1.16633  0.06556 2.347   0.0189 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                    exp(coef) exp(-coef) lower .95 upper .95
exposure_start_day1     1.289     0.7757    1.0542     1.576
exposure_start_day2     1.117     0.8951    0.9539     1.309
exposure_start_day3     1.154     0.8665    1.0111     1.317
exposure_start_day4     1.113     0.8984    0.9734     1.273
exposure_start_day5     1.166     0.8574    1.0257     1.326

Concordance= 0.725  (se = 0.008 )
Likelihood ratio test= 12.52  on 5 df,   p=0.03
Wald test            = 12.65  on 5 df,   p=0.03
Score (logrank) test = 12.69  on 5 df,   p=0.03
 

7.1.2 SCCS nevakcinovaných

I v případě neočkovaných jsou rizika po celou dobu relativně vyšší. Díky většímu množství dat je také možné zhodnotit tyto výsledky na vyšší hladině statistické významnosti.

combined_unvaccinated_prescribed <- dbGetQuery(con,
                                               "
                                    SELECT
                                    p.client_id,
                                    DATE '2021-06-15' AS vaccination_date,
                                    p.date AS event_date,
                                    'CPZP' as ins_company
                                    FROM cpzp_prescriptions p
                                    LEFT JOIN cpzp_clients c ON P.client_id = c.client_id 
                                    LEFT JOIN cpzp_ins_start_end s ON c.client_id = s.client_id
                                    WHERE 1 = 1
                                    AND s.ins_start_date < DATE '2021-06-15'- INTERVAL '1 year'
                                    AND s.ins_end_date > DATE '2021-06-15' + INTERVAL '1 year'
                                    AND 2021 - c.birthyear BETWEEN 16 AND 39
                                    AND p.date BETWEEN DATE '2021-06-15' - INTERVAL '1 year' AND DATE '2021-06-15' + INTERVAL '1 year'
                                    AND p.presc_order = 1
                                    AND c.n_vacc = 0
                                    
                                    UNION
                                    
                                    SELECT
                                    p.client_id,
                                    DATE '2021-06-15' AS vaccination_date,
                                    p.date AS event_date,
                                    'OZP' as ins_company
                                    FROM ozp_prescriptions p
                                    LEFT JOIN ozp_clients c ON P.client_id = c.client_id 
                                    LEFT JOIN ozp_ins_start_end s ON c.client_id = s.client_id
                                    WHERE 1 = 1
                                    AND s.ins_start_date < DATE '2021-06-15' - INTERVAL '1 year'
                                    AND s.ins_end_date > DATE '2021-06-15' + INTERVAL '1 year'
                                    AND 2021 - c.birthyear BETWEEN 16 AND 39
                                    AND p.date BETWEEN DATE '2021-06-15' - INTERVAL '1 year' AND DATE '2021-06-15' + INTERVAL '1 year'
                                    AND p.presc_order = 1
                                    AND c.n_vacc = 0

                                               ") %>%
  mutate(client_id = row_number()) %>%
  transmute(
    person_id = client_id,
    exposure_start_date = vaccination_date,
    event_date,
    observation_start_date = vaccination_date - days(365),
    observation_end_date = vaccination_date + days(365)
  )

sccs_unvaccinated <- combined_unvaccinated_prescribed %>%
  mutate(
    event_day = as.integer(event_date - observation_start_date),
    exposure_start_day = as.integer(exposure_start_date - observation_start_date),
    observation_start_day = 0,
    observation_end_day = as.integer(observation_end_date - observation_start_date)
  )

sccs_result_unvaccinated <- standardsccs(
  formula = event_day ~ exposure_start_day,
  indiv = person_id,
  astart = observation_start_day,
  aend = observation_end_day,
  aevent = event_day,
  adrug = exposure_start_day,
  aedrug = exposure_start_day + 365,
  expogrp = list(c(0, 30, 90, 180, 270)),
  data = sccs_unvaccinated
)

cat("```\n", paste(capture.output(print(sccs_result_unvaccinated)), collapse = "\n"), "\n```")
 Call:
coxph(formula = Surv(rep(1, 39084L), event) ~ exposure_start_day + 
    strata(indivL) + offset(log(interval)), data = chopdat, method = "exact")

  n= 39084, number of events= 6514 

                       coef exp(coef) se(coef)     z Pr(>|z|)    
exposure_start_day1 0.23835   1.26916  0.05838 4.083 4.45e-05 ***
exposure_start_day2 0.08343   1.08701  0.04608 1.810  0.07024 .  
exposure_start_day3 0.08523   1.08897  0.03900 2.185  0.02886 *  
exposure_start_day4 0.01957   1.01976  0.04003 0.489  0.62489    
exposure_start_day5 0.11229   1.11884  0.03763 2.984  0.00285 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                    exp(coef) exp(-coef) lower .95 upper .95
exposure_start_day1     1.269     0.7879    1.1319     1.423
exposure_start_day2     1.087     0.9200    0.9931     1.190
exposure_start_day3     1.089     0.9183    1.0088     1.175
exposure_start_day4     1.020     0.9806    0.9428     1.103
exposure_start_day5     1.119     0.8938    1.0393     1.204

Concordance= 0.731  (se = 0.005 )
Likelihood ratio test= 24.44  on 5 df,   p=2e-04
Wald test            = 25.16  on 5 df,   p=1e-04
Score (logrank) test = 25.22  on 5 df,   p=1e-04
 

7.1.3 Porovnání rizik

Průběh rizik je pro obě skupiny v zásadě srovnatelný s vyšším stanoveným rizikem pro očkovanou populaci. Ani jeden z rozdílů však není statisticky signifikantní.

coef_unvax <- sccs_result_unvaccinated$coefficients
conf_unvax <- sccs_result_unvaccinated$conf.int

df_unvax <- data.frame(
  term = rownames(coef_unvax),
  logHR = coef_unvax[, "coef"],
  SE = coef_unvax[, "se(coef)"],
  HR = coef_unvax[, "exp(coef)"],
  conf.low = conf_unvax[, "lower .95"],
  conf.high = conf_unvax[, "upper .95"]
)

coef_vax <- sccs_result_vaccinated$coefficients
conf_vax <- sccs_result_vaccinated$conf.int

df_vax <- data.frame(
  term = rownames(coef_vax),
  logHR = coef_vax[, "coef"],
  SE = coef_vax[, "se(coef)"],
  HR = coef_vax[, "exp(coef)"],
  conf.low = conf_vax[, "lower .95"],
  conf.high = conf_vax[, "upper .95"]
)

df_unvax <- df_unvax %>% mutate(group = "Bez Vakcinace")
df_vax <- df_vax %>% mutate(group = "Vakcinace")

combined <- bind_rows(df_vax, df_unvax)

combined$period <- as.integer(gsub("exposure_start_day", "", combined$term))

start_days <- c(0, 30, 90, 180, 270)
end_days   <- c(29, 89, 179, 269, Inf)

period_labels <- ifelse(
  is.infinite(end_days),
  paste0(start_days, "+"),
  paste0(start_days, "-", end_days)
)

combined$period_label <- factor(combined$period,
                                levels = seq_along(period_labels),
                                labels = period_labels)

ggplot(combined, aes(x = period_label, y = HR, color = group)) +
  geom_point(position = position_dodge(width = 0.6), size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2,
                position = position_dodge(width = 0.6)) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  labs(
    x = "Období expozice (dny)",
    y = "Poměr rizik (HR)",
    title = "Porovnání SCCS: očkovaní vs neočkovaní"
  ) +
  theme_minimal()

vax_df <- combined %>%
  filter(group == "Vakcinace") %>%
  rename_with(~ paste0(., "_vax"), -c(period, period_label))

unvax_df <- combined %>%
  filter(group == "Bez Vakcinace") %>%
  rename_with(~ paste0(., "_unvax"), -c(period, period_label))

df_diff <- inner_join(vax_df, unvax_df, by = c("period", "period_label"))

df_diff <- df_diff %>%
  mutate(
    logHR_diff = logHR_vax - logHR_unvax,
    se_diff = sqrt(SE_vax^2 + SE_unvax^2),
    ci_low = logHR_diff - 1.96 * se_diff,
    ci_high = logHR_diff + 1.96 * se_diff
  )

ggplot(df_diff, aes(x = period_label, y = logHR_diff)) +
  geom_point(color = "blue", size = 3) +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.2, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "Rozdíl v log(HR) (Očkovaní - Neočkovaní)",
    x = "Období expozice",
    y = "Rozdíl v log(HR)"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

common_terms <- intersect(df_vax$term, df_unvax$term)

did_results <- data.frame(
  term = character(),
  period_label = character(),
  logHR_vax = numeric(),
  logHR_unvax = numeric(),
  DiD_logHR = numeric(),
  DiD_SE = numeric(),
  Z = numeric(),
  P_value = numeric(),
  stringsAsFactors = FALSE
)

for (term in common_terms) {
  logHR_vax <- df_vax %>% filter(term == !!term) %>% pull(logHR)
  SE_vax    <- df_vax %>% filter(term == !!term) %>% pull(SE)
  
  logHR_unvax <- df_unvax %>% filter(term == !!term) %>% pull(logHR)
  SE_unvax    <- df_unvax %>% filter(term == !!term) %>% pull(SE)
  
  DiD_logHR <- logHR_vax - logHR_unvax
  DiD_SE <- sqrt(SE_vax^2 + SE_unvax^2)
  Z <- DiD_logHR / DiD_SE
  P_value <- 2 * pnorm(-abs(Z))
  
  period_idx <- as.integer(gsub("exposure_start_day", "", term))
  period_label <- period_labels[period_idx]
  
  did_results <- rbind(did_results, data.frame(
    term = term,
    period_label = period_label,
    logHR_vax = logHR_vax,
    logHR_unvax = logHR_unvax,
    DiD_logHR = DiD_logHR,
    DiD_SE = DiD_SE,
    Z = Z,
    P_value = P_value
  ))
}

did_results %>%
  transmute(Období = period_label, `log(HR) Vakcinace` = logHR_vax, `log(HR) Bez vakcinace` = logHR_unvax, `DiD log(HR)` = DiD_logHR, p_value = P_value) %>%
  kable(format = "html", caption = "Porovnání rizik očkovaných a neočkovaných (SCCS)") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání rizik očkovaných a neočkovaných (SCCS)
Období log(HR) Vakcinace log(HR) Bez vakcinace DiD log(HR) p_value
0-29 0.2539458 0.2383525 0.0155933 0.8949331
30-89 0.1108449 0.0834299 0.0274150 0.7678754
90-179 0.1432802 0.0852301 0.0580501 0.4564462
180-269 0.1071752 0.0195718 0.0876034 0.2692642
270+ 0.1538623 0.1122887 0.0415736 0.5823525

7.2 Předpisy celkem

7.2.1 SCCS vakcinovaných

Vyšší množství dat pro úhrn předpisů celkem umožňuje čistší pohled na zhodnocená rizika. Ta jsou stejně jako u prvopředpisů po celé sledované období vyšší než baseline.

combined_vaccinated_prescribed <- dbGetQuery(con,
                                    "
                                    SELECT
                                    v.client_id,
                                    v.date AS vaccination_date,
                                    p.date AS event_date,
                                    d.drug_form,
                                    'CPZP' as ins_company
                                    FROM cpzp_vaccinations v
                                    LEFT JOIN cpzp_clients c ON v.client_id = c.client_id
                                    LEFT JOIN cpzp_ins_start_end s ON c.client_id = s.client_id
                                    INNER JOIN cpzp_prescriptions p ON c.client_id = p.client_id
                                                                    AND p.date BETWEEN v.date - INTERVAL '1 year' AND v.date + INTERVAL '1 year'
                                    LEFT JOIN cpzp_drugs_catalog d ON p.presc_id  = d.presc_id 
                                    WHERE v.event_order = 1
                                    AND s.ins_start_date < v.date - INTERVAL '1 year'
                                    AND s.ins_end_date > v.date + INTERVAL '1 year'
                                    AND v.date BETWEEN DATE '2021-06-15' - INTERVAL '7 days'
                                                       AND DATE '2021-06-15' + INTERVAL '7 days'
                                    AND 2021 - c.birthyear BETWEEN 16 AND 39
                                    
                                    UNION
                                    
                                    SELECT
                                    v.client_id,
                                    v.date AS vaccination_date,
                                    p.date AS event_date,
                                    d.drug_form,
                                    'OZP' AS ins_company
                                    FROM OZP_vaccinations v
                                    LEFT JOIN ozp_clients c ON v.client_id = c.client_id
                                    LEFT JOIN ozp_ins_start_end s ON c.client_id = s.client_id
                                    INNER JOIN ozp_prescriptions p ON c.client_id = p.client_id
                                                                    AND p.date BETWEEN v.date - INTERVAL '1 year' AND v.date + INTERVAL '1 year'
                                    LEFT JOIN ozp_drugs_catalog d ON p.presc_id  = d.presc_id 
                                    WHERE v.event_order = 1
                                    AND s.ins_start_date < v.date - INTERVAL '1 year'
                                    AND s.ins_end_date > v.date + INTERVAL '1 year'
                                    AND v.date BETWEEN DATE '2021-06-15' - INTERVAL '7 days'
                                                       AND DATE '2021-06-15' + INTERVAL '7 days'
                                    AND 2021 - c.birthyear BETWEEN 16 AND 39
                                    ") %>%
  mutate(client_id = as.numeric(factor(client_id)),
         drug_form = tolower(drug_form),
         inj_or_inf = str_detect(drug_form, "injekční|infuzní"))

injectable_infusional <- combined_vaccinated_prescribed %>%
  filter(inj_or_inf) %>%
  arrange(client_id, event_date) %>%
  group_by(client_id) %>%
  mutate(
    day_diff = as.numeric(difftime(event_date, lag(event_date, default = first(event_date)), units = "days")),
    new_window = is.na(day_diff) | day_diff > 10,
    window_id = cumsum(new_window)
  ) %>%
  group_by(client_id, window_id) %>%
  slice(1) %>%
  ungroup()

other_drugs <- combined_vaccinated_prescribed %>%
  filter(!inj_or_inf)

combined_vaccinated_prescribed <- bind_rows(injectable_infusional, other_drugs) %>%
  arrange(client_id, event_date) %>%
  transmute(
    person_id = client_id,
    exposure_start_date = vaccination_date,
    event_date,
    observation_start_date = vaccination_date - days(365),
    observation_end_date = vaccination_date + days(365)
  )

sccs_vaccinated <- combined_vaccinated_prescribed %>%
  mutate(
    event_day = as.integer(event_date - observation_start_date),
    exposure_start_day = as.integer(exposure_start_date - observation_start_date),
    observation_start_day = 0,
    observation_end_day = as.integer(observation_end_date - observation_start_date)
  )

sccs_result_vaccinated <- standardsccs(
  formula = event_day ~ exposure_start_day,
  indiv = person_id,
  astart = observation_start_day,
  aend = observation_end_day,
  aevent = event_day,
  adrug = exposure_start_day,
  aedrug = exposure_start_day + 365,
  expogrp = list(c(0, 30, 90, 180, 270)),
  data = sccs_vaccinated
)

cat("```\n", paste(capture.output(print(sccs_result_vaccinated)), collapse = "\n"), "\n```")
 Call:
coxph(formula = Surv(rep(1, 33414L), event) ~ exposure_start_day + 
    strata(indivL) + offset(log(interval)), data = chopdat, method = "exact")

  n= 33414, number of events= 5569 

                       coef exp(coef) se(coef)     z Pr(>|z|)    
exposure_start_day1 0.09611   1.10088  0.06856 1.402  0.16094    
exposure_start_day2 0.11954   1.12698  0.04995 2.393  0.01671 *  
exposure_start_day3 0.11603   1.12303  0.04241 2.736  0.00622 ** 
exposure_start_day4 0.10327   1.10879  0.04262 2.423  0.01539 *  
exposure_start_day5 0.27407   1.31531  0.03895 7.036 1.98e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                    exp(coef) exp(-coef) lower .95 upper .95
exposure_start_day1     1.101     0.9084    0.9625     1.259
exposure_start_day2     1.127     0.8873    1.0219     1.243
exposure_start_day3     1.123     0.8904    1.0335     1.220
exposure_start_day4     1.109     0.9019    1.0199     1.205
exposure_start_day5     1.315     0.7603    1.2186     1.420

Concordance= 0.732  (se = 0.005 )
Likelihood ratio test= 50.28  on 5 df,   p=1e-09
Wald test            = 51.74  on 5 df,   p=6e-10
Score (logrank) test = 51.97  on 5 df,   p=5e-10
 

7.2.2 SCCS nevakcinovaných

Je patrné že průběh rizik pro neočkované téměř kopíruje průběh rizik očkovaných jedinců.

combined_unvaccinated_prescribed <- dbGetQuery(con,
                                               "
                                    SELECT
                                    p.client_id,
                                    DATE '2021-06-15' AS vaccination_date,
                                    p.date AS event_date,
                                    d.drug_form,
                                    'CPZP' as ins_company
                                    FROM cpzp_prescriptions p
                                    LEFT JOIN cpzp_clients c ON P.client_id = c.client_id 
                                    LEFT JOIN cpzp_ins_start_end s ON c.client_id = s.client_id
                                    LEFT JOIN cpzp_drugs_catalog d ON p.presc_id  = d.presc_id 
                                    WHERE 1 = 1
                                    AND s.ins_start_date < DATE '2021-06-15' - INTERVAL '1 year'
                                    AND s.ins_end_date > DATE '2021-06-15' + INTERVAL '1 year'
                                    AND 2021 - c.birthyear BETWEEN 16 AND 39
                                    AND p.date BETWEEN DATE '2021-06-15' - INTERVAL '1 year' AND DATE '2021-06-15' + INTERVAL '1 year'
                                    AND c.n_vacc = 0
                                    
                                    UNION
                                    
                                    SELECT
                                    p.client_id,
                                    DATE '2021-06-15' AS vaccination_date,
                                    p.date AS event_date,
                                    d.drug_form,
                                    'OZP' as ins_company
                                    FROM ozp_prescriptions p
                                    LEFT JOIN ozp_clients c ON P.client_id = c.client_id 
                                    LEFT JOIN ozp_ins_start_end s ON c.client_id = s.client_id
                                    LEFT JOIN ozp_drugs_catalog d ON p.presc_id  = d.presc_id 
                                    WHERE 1 = 1
                                    AND s.ins_start_date < DATE '2021-06-15' - INTERVAL '1 year'
                                    AND s.ins_end_date > DATE '2021-06-15' + INTERVAL '1 year'
                                    AND 2021 - c.birthyear BETWEEN 16 AND 39
                                    AND p.date BETWEEN DATE '2021-06-15' - INTERVAL '1 year' AND DATE '2021-06-15' + INTERVAL '1 year'
                                    AND c.n_vacc = 0

                                               ") %>%
  mutate(client_id = as.numeric(factor(client_id)),
         drug_form = tolower(drug_form),
         inj_or_inf = str_detect(drug_form, "injekční|infuzní"))

injectable_infusional <- combined_unvaccinated_prescribed %>%
  filter(inj_or_inf) %>%
  arrange(client_id, event_date) %>%
  group_by(client_id) %>%
  mutate(
    day_diff = as.numeric(difftime(event_date, lag(event_date, default = first(event_date)), units = "days")),
    new_window = is.na(day_diff) | day_diff > 10,
    window_id = cumsum(new_window)
  ) %>%
  group_by(client_id, window_id) %>%
  slice(1) %>%
  ungroup()

other_drugs <- combined_unvaccinated_prescribed %>%
  filter(!inj_or_inf)

combined_unvaccinated_prescribed <- bind_rows(injectable_infusional, other_drugs) %>%
  arrange(client_id, event_date) %>%
  transmute(
    person_id = client_id,
    exposure_start_date = vaccination_date,
    event_date,
    observation_start_date = vaccination_date - days(365),
    observation_end_date = vaccination_date + days(365)
  )

sccs_unvaccinated <- combined_unvaccinated_prescribed %>%
  mutate(
    event_day = as.integer(event_date - observation_start_date),
    exposure_start_day = as.integer(exposure_start_date - observation_start_date),
    observation_start_day = 0,
    observation_end_day = as.integer(observation_end_date - observation_start_date)
  )

sccs_result_unvaccinated <- standardsccs(
  formula = event_day ~ exposure_start_day,
  indiv = person_id,
  astart = observation_start_day,
  aend = observation_end_day,
  aevent = event_day,
  adrug = exposure_start_day,
  aedrug = exposure_start_day + 365,
  expogrp = list(c(0, 30, 90, 180, 270)),
  data = sccs_unvaccinated
)

cat("```\n", paste(capture.output(print(sccs_result_unvaccinated)), collapse = "\n"), "\n```")
 Call:
coxph(formula = Surv(rep(1, 109116L), event) ~ exposure_start_day + 
    strata(indivL) + offset(log(interval)), data = chopdat, method = "exact")

  n= 109116, number of events= 18186 

                       coef exp(coef) se(coef)     z Pr(>|z|)    
exposure_start_day1 0.19387   1.21394  0.03608 5.374 7.71e-08 ***
exposure_start_day2 0.05440   1.05591  0.02826 1.925 0.054178 .  
exposure_start_day3 0.14992   1.16174  0.02303 6.510 7.50e-11 ***
exposure_start_day4 0.08169   1.08512  0.02365 3.454 0.000552 ***
exposure_start_day5 0.19609   1.21663  0.02207 8.885  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                    exp(coef) exp(-coef) lower .95 upper .95
exposure_start_day1     1.214     0.8238     1.131     1.303
exposure_start_day2     1.056     0.9471     0.999     1.116
exposure_start_day3     1.162     0.8608     1.110     1.215
exposure_start_day4     1.085     0.9216     1.036     1.137
exposure_start_day5     1.217     0.8219     1.165     1.270

Concordance= 0.731  (se = 0.003 )
Likelihood ratio test= 112  on 5 df,   p=<2e-16
Wald test            = 113.6  on 5 df,   p=<2e-16
Score (logrank) test = 113.8  on 5 df,   p=<2e-16
 

7.2.3 Porovnání rizik

Porovnatelnost mezi skupinami je téměř nulová. Jediný rozdíl je v období 270-365, kde je vyšší riziko pro očkované s p-hodnotou Waldova testu těsně nad úrovní 0,08.

coef_unvax <- sccs_result_unvaccinated$coefficients
conf_unvax <- sccs_result_unvaccinated$conf.int

df_unvax <- data.frame(
  term = rownames(coef_unvax),
  logHR = coef_unvax[, "coef"],
  SE = coef_unvax[, "se(coef)"],
  HR = coef_unvax[, "exp(coef)"],
  conf.low = conf_unvax[, "lower .95"],
  conf.high = conf_unvax[, "upper .95"]
)

coef_vax <- sccs_result_vaccinated$coefficients
conf_vax <- sccs_result_vaccinated$conf.int

df_vax <- data.frame(
  term = rownames(coef_vax),
  logHR = coef_vax[, "coef"],
  SE = coef_vax[, "se(coef)"],
  HR = coef_vax[, "exp(coef)"],
  conf.low = conf_vax[, "lower .95"],
  conf.high = conf_vax[, "upper .95"]
)

df_unvax <- df_unvax %>% mutate(group = "Unvaccinated")
df_vax <- df_vax %>% mutate(group = "Vaccinated")

combined <- bind_rows(df_vax, df_unvax)

combined$period <- as.integer(gsub("exposure_start_day", "", combined$term))

start_days <- c(0, 30, 90, 180, 270)
end_days   <- c(29, 89, 179, 269, Inf)

period_labels <- ifelse(
  is.infinite(end_days),
  paste0(start_days, "+"),
  paste0(start_days, "-", end_days)
)

combined$period_label <- factor(combined$period,
                                levels = seq_along(period_labels),
                                labels = period_labels)

ggplot(combined, aes(x = period_label, y = HR, color = group)) +
  geom_point(position = position_dodge(width = 0.6), size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2,
                position = position_dodge(width = 0.6)) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  labs(
    x = "Období expozice (dny)",
    y = "Poměr rizik (HR)",
    title = "Porovnání SCCS: očkovaní vs neočkovaní"
  ) +
  theme_minimal()

vax_df <- combined %>%
  filter(group == "Vaccinated") %>%
  rename_with(~ paste0(., "_vax"), -c(period, period_label))

unvax_df <- combined %>%
  filter(group == "Unvaccinated") %>%
  rename_with(~ paste0(., "_unvax"), -c(period, period_label))

df_diff <- inner_join(vax_df, unvax_df, by = c("period", "period_label"))

df_diff <- df_diff %>%
  mutate(
    logHR_diff = logHR_vax - logHR_unvax,
    se_diff = sqrt(SE_vax^2 + SE_unvax^2),
    ci_low = logHR_diff - 1.96 * se_diff,
    ci_high = logHR_diff + 1.96 * se_diff
  )

ggplot(df_diff, aes(x = period_label, y = logHR_diff)) +
  geom_point(color = "blue", size = 3) +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.2, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "Rozdíl v log(HR) (Očkovaní - Neočkovaní)",
    x = "Období expozice",
    y = "Rozdíl v log(HR)"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

common_terms <- intersect(df_vax$term, df_unvax$term)

did_results <- data.frame(
  term = character(),
  period_label = character(),
  logHR_vax = numeric(),
  logHR_unvax = numeric(),
  DiD_logHR = numeric(),
  DiD_SE = numeric(),
  Z = numeric(),
  P_value = numeric(),
  stringsAsFactors = FALSE
)

for (term in common_terms) {
  logHR_vax <- df_vax %>% filter(term == !!term) %>% pull(logHR)
  SE_vax    <- df_vax %>% filter(term == !!term) %>% pull(SE)
  
  logHR_unvax <- df_unvax %>% filter(term == !!term) %>% pull(logHR)
  SE_unvax    <- df_unvax %>% filter(term == !!term) %>% pull(SE)
  
  DiD_logHR <- logHR_vax - logHR_unvax
  DiD_SE <- sqrt(SE_vax^2 + SE_unvax^2)
  Z <- DiD_logHR / DiD_SE
  P_value <- 2 * pnorm(-abs(Z))
  
  period_idx <- as.integer(gsub("exposure_start_day", "", term))
  period_label <- period_labels[period_idx]
  
  did_results <- rbind(did_results, data.frame(
    term = term,
    period_label = period_label,
    logHR_vax = logHR_vax,
    logHR_unvax = logHR_unvax,
    DiD_logHR = DiD_logHR,
    DiD_SE = DiD_SE,
    Z = Z,
    P_value = P_value
  ))
}

did_results %>%
  transmute(Období = period_label, `log(HR) Vakcinace` = logHR_vax, `log(HR) Bez vakcinace` = logHR_unvax, `DiD log(HR)` = DiD_logHR, p_value = P_value) %>%
  kable(format = "html", caption = "Porovnání rizik (SCCS) očkovaných a neočkovaných") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání rizik (SCCS) očkovaných a neočkovaných
Období log(HR) Vakcinace log(HR) Bez vakcinace DiD log(HR) p_value
0-29 0.0961135 0.1938689 -0.0977554 0.2070137
30-89 0.1195438 0.0544028 0.0651410 0.2563519
90-179 0.1160288 0.1499191 -0.0338902 0.4825014
180-269 0.1032718 0.0816914 0.0215803 0.6579593
270+ 0.2740703 0.1960854 0.0779850 0.0815362

7.3 Shrnutí

Krátkodobá rizika jsou ve všech případech vyšší oproti baseline a v zásadě srovnatelná. Jediný částečně významný rozdíl je u celkového úhrnu předpisů v období 270-365 s vyšším rizikem pro očkované.

8 Věková kategorie 40-59

8.1 Dlouhodobý efekt

8.1.1 Distribuce věku a pohlaví

unique_clients <- prescriptions %>%
  filter(!is.na(age_2021), !is.na(sex), !is.na(vaccinated)) %>%
  distinct(client_id, .keep_all = TRUE)  

age_data <- unique_clients %>%
  group_by(vaccinated, age_2021) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(vaccinated) %>%
  mutate(prop = count / sum(count),
         variable = "Age",
         value = as.character(age_2021)) %>%
  ungroup()

sex_data <- unique_clients %>%
  group_by(vaccinated, sex) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(vaccinated) %>%
  mutate(prop = count / sum(count),
         variable = "Sex",
         value = sex) %>%
  ungroup()

Stáří populace a poměr dle pohlaví jsou v zásadě srovnatelné.

ggplot(unique_clients, aes(x = age_2021, fill = factor(vaccinated))) +
  geom_histogram(aes(y = after_stat(density)), position = "identity", alpha = 0.4, binwidth = 1) +
  scale_fill_brewer(palette = "Set1", name = "Vaccinated") +
  scale_x_continuous(breaks = seq(floor(min(unique_clients$age)), ceiling(max(unique_clients$age)), by = 1)) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    title = "Věkové rozložení podle statusu vakcinace",
    x = "Věk", y = "%",
    fill = "Vakcinace"
  ) +
  theme_minimal()

ggplot(sex_data, aes(x = sex, y = prop, color = vaccinated, group = vaccinated)) +
  geom_line(size = 1) +
  geom_point(size = 2) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(
    title = "Pohlaví podle vakcinace",
    x = "Pohlaví",
    y = "%",
    color = "Vakcinace"
  ) +
  theme_minimal()

8.1.2 Předpověď vs realita

prescriptions_clean <- prescriptions %>%
  filter(!is.na(presc_date)) %>%
  mutate(presc_date = as.Date(presc_date)) %>%
  arrange(presc_date) %>%
  filter(presc_date >= as.Date("2015-01-01") & presc_date <= as.Date("2023-12-31")) %>%
  mutate(count = 1)  


results_vacc <- prescriptions_clean %>%
  filter(!is.na(vaccinated)) %>%
  group_by(vaccinated) %>%
  group_split() %>%
  map(~ {
    group_val <- unique(.x$vaccinated)
    analyze_group(df = filter(prescriptions_clean, vaccinated == group_val), 
                  group_id = group_val, 
                  group_type = "Vaccination")
  }) %>%
  compact() %>%
  bind_rows()

results_sex <- prescriptions_clean %>%
  filter(!is.na(sex)) %>%
  group_by(sex) %>%
  group_split() %>%
  map(~ {
    group_val <- unique(.x$sex)
    analyze_group(df = filter(prescriptions_clean, sex == group_val), 
                  group_id = group_val, 
                  group_type = "Sex")
  }) %>%
  compact() %>%
  bind_rows()

results_all <- bind_rows(results_vacc, results_sex) %>%
  arrange(p_value)

plot_data_all <- results_all %>%
  unnest(plot_data, names_sep = "_")

ETS model vykazuje odchylky napříč všemi skupinami s nejvýraznějším rozdílem u neočkované populace.

results_all %>%
  mutate(Skupina = group, RMSE = rmse, rRMSE = rrmse) %>%
  select(
    Skupina, RMSE, rRMSE, t_stat, p_value
  ) %>%
  kable(format = "html", caption = "Porovnání předpisovosti v čase (ETS)") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání předpisovosti v čase (ETS)
Skupina RMSE rRMSE t_stat p_value
Bez vakcinace 245.4343 0.1401934 -5.0840885 0.0000063
M 336.2235 0.1120938 -2.8601803 0.0062977
Vakcinace 704.0997 0.1028859 -0.7368828 0.4648563
F 576.8500 0.1031065 0.0795498 0.9369331
ggplot(plot_data_all, aes(x = plot_data_month)) +
  geom_line(aes(y = plot_data_prescriptions), color = "blue") +
  geom_line(aes(y = plot_data_predicted), color = "red", linetype = "dashed") +
  geom_vline(xintercept = as.Date("2020-01-01"), linetype = "dotted") +
  facet_wrap(~ paste(plot_data_group_type, plot_data_group, sep = ": "), scales = "free_y") +
  labs(
    title = "Realita vs Předpověď (ETS)",
    subtitle = "Modrá = Realita, Červená = Předpověď",
    x = "Datum", y = "Počet předpisů"
  ) +
  theme_minimal()

8.1.3 Prvopředpisovost

first_presc <- dbGetQuery(con,
                          "
                          SELECT
                            c.client_id,
                            c.n_vacc,
                            c.n_presc,
                            c.sex,
                            c.birthyear,
                            FIRST(p.date) OVER (PARTITION BY c.client_id ORDER BY p.date ASC) AS first_presc_date,
                            'CPZP' AS ins_company
                          FROM
                            cpzp_clients c
                          LEFT JOIN cpzp_ins_start_end i ON c.client_id = i.client_id
                          LEFT JOIN cpzp_prescriptions p ON c.client_id = p.client_id AND p.presc_order = 1
                          WHERE
                            i.ins_start_date < '2015-01-01' AND i.ins_end_date = '9999-12-31'
                            AND 2021 - c.birthyear BETWEEN 40 AND 59
                            
                          UNION
                          
                          SELECT
                            c.client_id,
                            c.n_vacc,
                            c.n_presc,
                            c.sex,
                            c.birthyear,
                            FIRST(p.date) OVER (PARTITION BY c.client_id ORDER BY p.date ASC) AS first_presc_date,
                            'OZP' AS ins_company
                          FROM
                            ozp_clients c
                          LEFT JOIN ozp_ins_start_end i ON c.client_id = i.client_id
                          LEFT JOIN ozp_prescriptions p ON c.client_id = p.client_id AND p.presc_order = 1
                          WHERE
                            i.ins_start_date < '2015-01-01' AND i.ins_end_date IS NULL
                            AND 2021 - c.birthyear BETWEEN 40 AND 59
                          ") %>%
  mutate(first_presc_yearmonth = floor_date(first_presc_date, unit = "month"),
         sex = ifelse(sex == "Z", "F", sex),
         vacc_status = ifelse(n_vacc > 0, "Vakcinace", "Bez vakcinace"),
         age_2021 = 2021 - birthyear)

monthly_first_presc <- first_presc %>%
  count(first_presc_yearmonth, name = "new_prescriptions") %>%
  arrange(first_presc_yearmonth) %>%
  mutate(
    cum_new_presc = cumsum(new_prescriptions),  
    total_population = nrow(first_presc),
    at_risk_population = total_population - lag(cum_new_presc, default = 0),
    rate = new_prescriptions / at_risk_population
  ) %>%
  mutate(
    period_group = case_when(
      first_presc_yearmonth < as.Date("2020-01-01") ~ "Pre-COVID",
      first_presc_yearmonth > as.Date("2021-12-31") ~ "Post-COVID",
      TRUE ~ "COVID-Shock"
    )
  )

monthly_clean <- monthly_first_presc %>%
  filter(first_presc_yearmonth >= as.Date("2018-01-01")) %>%
  filter(period_group %in% c("Pre-COVID", "Post-COVID"))

ttest <- t.test(rate ~ period_group, data = monthly_clean)

Průměrná prvopředpisovost ve srovnání po vs před je mírně nižší. Samotný rozdíl průměrů však není statisticky významný.

ggplot(monthly_first_presc %>% filter(first_presc_yearmonth >= as.Date("2018-01-01")),
       aes(x = first_presc_yearmonth, y = rate * 1000)) +
  geom_line(color = "steelblue", size = 1) +
  annotate("rect",
           xmin = as.Date("2020-01-01"), xmax = as.Date("2021-12-31"),
           ymin = -Inf, ymax = Inf, alpha = 0.2, fill = "gray50") +
  scale_y_continuous(name = "Počet prvopředpisů na 1000 osob bez předpisu", limits = c(0, 10)) +
  labs(
    title = "Prvopředpisovost v čase",
    subtitle = "Šedá = covidové období (vyřazeno z t-testu)",
    x = "Datum"
  ) +
  theme_minimal()

cat("```\n", paste(capture.output(ttest), collapse = "\n"), "\n```")
 
    Welch Two Sample t-test

data:  rate by period_group
t = -1.1168, df = 45.53, p-value = 0.27
alternative hypothesis: true difference in means between group Post-COVID and group Pre-COVID is not equal to 0
95 percent confidence interval:
 -0.0004915077  0.0001407987
sample estimates:
mean in group Post-COVID  mean in group Pre-COVID 
             0.005124781              0.005300136 
 
total_by_vacc <- first_presc %>%
  count(vacc_status, name = "total_population")

monthly_by_vacc <- first_presc %>%
  count(vacc_status, first_presc_yearmonth, name = "new_prescriptions") %>%
  arrange(vacc_status, first_presc_yearmonth) %>%
  group_by(vacc_status) %>%
  mutate(
    cum_prescriptions = cumsum(new_prescriptions),
    total_population = total_by_vacc$total_population[match(vacc_status, total_by_vacc$vacc_status)],
    at_risk_population = total_population - lag(cum_prescriptions, default = 0),
    rate = new_prescriptions / at_risk_population
  ) %>%
  ungroup() %>%
  mutate(
    period_group = case_when(
      first_presc_yearmonth < as.Date("2020-01-01") ~ "Pre-COVID",
      first_presc_yearmonth > as.Date("2021-12-31") ~ "Post-COVID",
      TRUE ~ "COVID-Shock"
    )
  )
  
monthly_by_vacc_clean <- monthly_by_vacc %>%
  filter(first_presc_yearmonth >= as.Date("2018-01-01")) %>%
  filter(period_group %in% c("Pre-COVID", "Post-COVID"))

first_presc_table <- monthly_by_vacc_clean %>%
  group_by(vacc_status) %>%
  summarise(
    p_value = t.test(rate ~ period_group)$p.value,
    mean_before = mean(rate[period_group == "Pre-COVID"]),
    mean_after = mean(rate[period_group == "Post-COVID"]),
    rate = mean_after / mean_before
  )

Obě skupiny vykazují mírné snížení v prvopředpisovosti. Rozdíl mezi skupinami je nevýznamný.

first_presc_table %>%
  transmute(Vakcinace = vacc_status, p_value, `Průměr před` = mean_before * 1000, `Průměr po` = mean_after * 1000, `Poměr` = rate) %>%
  kable(format = "html", caption = "Porovnání prvopředpisovosti v čase podle vakcinace") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání prvopředpisovosti v čase podle vakcinace
Vakcinace p_value Průměr před Průměr po Poměr
Bez vakcinace 0.0606129 4.325463 4.072212 0.9414511
Vakcinace 0.4742458 5.642638 5.518322 0.9779684
ggplot(monthly_by_vacc %>% filter(first_presc_yearmonth >= as.Date("2018-01-01")),
       aes(x = first_presc_yearmonth, y = rate * 1000, color = vacc_status)) +
  geom_line(size = 1) +
  annotate("rect",
           xmin = as.Date("2020-01-01"), xmax = as.Date("2021-12-31"),
           ymin = -Inf, ymax = Inf, alpha = 0.2, fill = "gray50") +
  scale_y_continuous(name = "Počet předpisů na 1000 osob bez předpisu", limits = c(0, 10)) +
  labs(
    title = "Prvopředpisovost v čase podle vakcinace",
    subtitle = "Šedá = covidové období (vyřazeno z t-testu)",
    x = "Datum",
    color = "Vakcinace"
  ) +
  theme_minimal()

monthly_rates <- monthly_by_vacc %>%
  filter(period_group %in% c("Pre-COVID", "Post-COVID")) %>%
  mutate(
    vaccinated = if_else(vacc_status == "Vakcinace", 1, 0),
    post_covid = if_else(first_presc_yearmonth >= as.Date("2022-01-01"), 1, 0),
    interaction = vaccinated * post_covid)

did_model <- lm(log(rate * 1000) ~ vaccinated * post_covid, data = monthly_rates)

cat("```\n", paste(capture.output(summary(did_model)), collapse = "\n"), "\n```")
 
Call:
lm(formula = log(rate * 1000) ~ vaccinated * post_covid, data = monthly_rates)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.35508 -0.10601 -0.02006  0.08872  0.75335 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)            1.58083    0.02292  68.980  < 2e-16 ***
vaccinated             0.24642    0.03241   7.603 2.12e-12 ***
post_covid            -0.18320    0.04287  -4.273 3.26e-05 ***
vaccinated:post_covid  0.05736    0.06063   0.946    0.345    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.1775 on 164 degrees of freedom
Multiple R-squared:  0.4203,    Adjusted R-squared:  0.4097 
F-statistic: 39.64 on 3 and 164 DF,  p-value: < 2.2e-16
 

8.1.4 Multipředpisovost

8.1.4.1 Míra retence (meziroční multipředpisovost)

presc_id_counts <- dbGetQuery(con,
                          "
                          SELECT
                            c.client_id,
                            c.n_vacc,
                            c.n_presc,
                            c.sex,
                            c.birthyear,
                            2021 - c.birthyear AS age_2021,
                            p.date as presc_date,
                            p.presc_order,
                            'CPZP' AS ins_company,
                            YEAR(FIRST(p.date) OVER (PARTITION BY c.client_id ORDER BY p.presc_order ASC)) AS first_presc_year
                          FROM
                            cpzp_clients c
                          LEFT JOIN cpzp_ins_start_end i ON c.client_id = i.client_id
                          LEFT JOIN cpzp_prescriptions p ON c.client_id = p.client_id
                          WHERE
                            i.ins_start_date < '2015-01-01' AND i.ins_end_date = '9999-12-31'
                            AND 2021 - c.birthyear BETWEEN 40 AND 59
                            
                          UNION
                          
                          SELECT
                            c.client_id,
                            c.n_vacc,
                            c.n_presc,
                            c.sex,
                            c.birthyear,
                            2021 - c.birthyear AS age_2021,
                            p.date as presc_date,
                            p.presc_order,
                            'OZP' AS ins_company,
                            YEAR(FIRST(p.date) OVER (PARTITION BY c.client_id ORDER BY p.presc_order ASC)) AS first_presc_year
                          FROM
                            ozp_clients c
                          LEFT JOIN ozp_ins_start_end i ON c.client_id = i.client_id
                          LEFT JOIN ozp_prescriptions p ON c.client_id = p.client_id
                          WHERE
                            i.ins_start_date < '2015-01-01' AND i.ins_end_date IS NULL
                            AND 2021 - c.birthyear BETWEEN 40 AND 59
                          ")%>%
  mutate(presc_year = as.integer(format(as.Date(presc_date), "%Y")),
         vacc_status = ifelse(n_vacc > 0, "Vakcinace", "Bez vakcinace"))

monthly_unique_clients <- presc_id_counts %>%
  mutate(yearmonth = floor_date(presc_date, unit = "month")) %>%
  distinct(client_id, yearmonth, vacc_status) %>%
  count(yearmonth, vacc_status, name = "n_unique_clients")

presc_years <- presc_id_counts %>%
  mutate(
    presc_year = year(presc_date)
  ) %>%
  distinct(client_id, first_presc_year, presc_year, vacc_status) %>%
  filter(first_presc_year >= 2017, first_presc_year <= 2023)


retention_data <- presc_years %>%
  mutate(year_offset = presc_year - first_presc_year) %>%
  filter(year_offset >= 0)


retention_summary <- retention_data %>%
  group_by(first_presc_year, year_offset, vacc_status) %>%
  summarise(n_clients = n_distinct(client_id), .groups = "drop")


cohort_sizes <- retention_data %>%
  filter(year_offset == 0) %>%
  group_by(first_presc_year, vacc_status) %>%
  summarise(cohort_size = n_distinct(client_id), .groups = "drop")


retention_summary <- retention_summary %>%
  left_join(cohort_sizes, by = c("first_presc_year", "vacc_status")) %>%
  mutate(retention_rate = n_clients / cohort_size)

retention_summary_filtered <- retention_summary %>%
  filter(year_offset > 0)

Míra meziročních předpisů stále potvrzuje trend z předchozích sledování. Je vyšší než u mladší populace a míra pro očkované je všeobecně vyšší.

ggplot(retention_summary_filtered, aes(x = year_offset, y = retention_rate, color = as.factor(first_presc_year))) +
  geom_line(size = 1.2) +
  geom_point(size = 2) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  scale_color_viridis_d(name = "Rok prvopředpisu") +
  labs(
    title = "Retence (meziroční předpisovost) podle vakcinace",
    x = "Počet let od prvopředpisu",
    y = "Míra multipředpisovosti"
  ) +
  facet_wrap(~ vacc_status) +
  theme_minimal(base_size = 13)

8.1.4.2 Roční multipředpisovost

multi_presc_yearly <- multi_presc %>%
  mutate(year = year(date)) %>%
  filter(year > 2017)


multi_clients_yearly <- multi_presc_yearly %>%
  group_by(year, vacc_status, client_id) %>%
  summarise(presc_count = n(), .groups = "drop") %>%
  mutate(is_multi = presc_count > 1)


multi_summary_yearly <- multi_clients_yearly %>%
  group_by(year, vacc_status) %>%
  summarise(
    total_clients = n(),
    multi_clients = sum(is_multi),
    multi_rate = (multi_clients / total_clients) * 1000,
    .groups = "drop"
  )


multi_summary_yearly <- multi_summary_yearly %>%
  mutate(
    period_group = case_when(
      year < 2020 ~ "Pre-COVID",
      year > 2021 ~ "Post-COVID",
      TRUE ~ "COVID-Shock"
    )
  )


multi_summary_clean <- multi_summary_yearly %>%
  filter(period_group %in% c("Pre-COVID", "Post-COVID"))


t_test_results <- multi_summary_clean %>%
  group_by(vacc_status) %>%
  summarise(
    p_value = t.test(multi_rate ~ period_group)$p.value,
    mean_before = mean(multi_rate[period_group == "Pre-COVID"]),
    mean_after = mean(multi_rate[period_group == "Post-COVID"]),
    rate = mean_after / mean_before,
    .groups = "drop"
  )

Roční multipředpisovost stoupá v obou kategoriích. Míra změny mezi skupinami je nevýznamná.

t_test_results %>%
  transmute(Vakcinace = vacc_status, p_value, `Průměr před` = mean_before, `Průměr po` = mean_after, `Poměr` = rate) %>%
  kable(format = "html", caption = "Porovnání roční multipředpisovosti v čase") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání roční multipředpisovosti v čase
Vakcinace p_value Průměr před Průměr po Poměr
Bez vakcinace 0.0892644 409.2620 434.4393 1.061519
Vakcinace 0.1074908 419.2872 438.8464 1.046649
ggplot(multi_summary_yearly, aes(x = year, y = multi_rate, color = vacc_status)) +
  geom_line(size = 0.6, alpha = 0.4) +
  stat_smooth(method = "loess", se = FALSE, size = 1.2) +
  annotate("rect",
           xmin = 2020, xmax = 2021,
           ymin = -Inf, ymax = Inf,
           alpha = 0.2, fill = "gray50") +
  scale_y_continuous(name = "Roční multipředpisovost na 1000 předepsaných klientů") +
  scale_x_continuous(breaks = unique(multi_summary_yearly$year)) +
  labs(
    title = "Roční multipředpisovost podle vakcinace",
    subtitle = "Šedá = covidové období (vyřazeno z t-testu)",
    x = "Rok",
    color = "Vakcinace"
  ) +
  theme_minimal()

did_data <- multi_summary_clean %>%
  mutate(
    vaccinated = if_else(vacc_status == "Vakcinace", 1, 0),
    post_covid = if_else(year >= 2022, 1, 0),
    interaction = vaccinated * post_covid
  )

did_model <- lm(log(multi_rate) ~ vaccinated * post_covid, data = did_data)
cat("```\n", paste(capture.output(summary(did_model)), collapse = "\n"), "\n```")
 
Call:
lm(formula = log(multi_rate) ~ vaccinated * post_covid, data = did_data)

Residuals:
        1         2         3         4         5         6         7         8 
-0.014567 -0.012360  0.014567  0.012360 -0.008741 -0.007212  0.008741  0.007212 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)            6.01425    0.01111 541.532 6.98e-11 ***
vaccinated             0.02423    0.01571   1.543    0.198    
post_covid             0.05977    0.01571   3.805    0.019 *  
vaccinated:post_covid -0.01412    0.02221  -0.636    0.559    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.01571 on 4 degrees of freedom
Multiple R-squared:  0.8636,    Adjusted R-squared:  0.7612 
F-statistic: 8.439 on 3 and 4 DF,  p-value: 0.03328
 

8.1.4.3 Měsíční multipředpisovost

Míra měsíční multipředpisovosti vykazuje zajímavé rozdělení trendů opačnými směry v období covidu, kdy všeobecné riziko bylo u vakcinovaných jedinců menší. Tento rozdíl je statisticky významný.

t_test_results %>%
  transmute(Vakcinace = vacc_status, p_value, `Průměr před` = mean_before, `Průměr po` = mean_after, `Poměr` = rate) %>%
  kable(format = "html", caption = "Porovnání měsíční multipředpisovosti v čase") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání měsíční multipředpisovosti v čase
Vakcinace p_value Průměr před Průměr po Poměr
Bez vakcinace 0.0002998 151.3781 162.6159 1.074237
Vakcinace 0.5848585 150.4276 151.5545 1.007491
ggplot(multi_summary_monthly, aes(x = year_month, y = multi_rate, color = vacc_status)) +
  geom_line(size = 0.5, alpha = 0.4) + 
  stat_smooth(method = "loess", se = FALSE, size = 1.2) +  
  annotate("rect",
           xmin = as.Date("2020-01-01"), xmax = as.Date("2021-12-31"),
           ymin = -Inf, ymax = Inf, alpha = 0.2, fill = "gray50") +
  scale_y_continuous(name = "Měsíční multipředpisovost na 1000 předepsaných klientů") +
  labs(
    title = "Měsíční multipředpisovost podle vakcinačního statusu",
    subtitle = "Šedá = covidové období (vyřazeno z t-testu)",
    x = "Datum",
    color = "Vakcinace"
  ) +
  theme_minimal()

did_data <- multi_summary_clean %>%
  mutate(
    vaccinated = if_else(vacc_status == "Vakcinace", 1, 0),
    post_covid = if_else(year_month >= as.Date("2022-01-01"), 1, 0),
    interaction = vaccinated * post_covid
  )

did_model <- lm(log(multi_rate) ~ vaccinated * post_covid, data = did_data)
cat("```\n", paste(capture.output(summary(did_model)), collapse = "\n"), "\n```")
 
Call:
lm(formula = log(multi_rate) ~ vaccinated * post_covid, data = did_data)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.122492 -0.038151 -0.003402  0.035210  0.183638 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)    
(Intercept)            5.01774    0.01129 444.295  < 2e-16 ***
vaccinated            -0.00548    0.01597  -0.343  0.73229    
post_covid             0.07192    0.01597   4.503 1.96e-05 ***
vaccinated:post_covid -0.06415    0.02259  -2.840  0.00555 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.05533 on 92 degrees of freedom
Multiple R-squared:  0.2555,    Adjusted R-squared:  0.2312 
F-statistic: 10.52 on 3 and 92 DF,  p-value: 5.127e-06
 

8.1.5 Způsob podání

Trendy nastavené v předchozích kohortách jsou v zásadě nezměněné i ve starší populaci.

prescriptions <- prescriptions %>%
  mutate(
    drug_form_std = case_when(
      is.na(drug_form) ~ NA_character_,
      str_detect(drug_form, "tableta|Tableta") ~ "Tableta",
      str_detect(drug_form, "tobolka|Tobolka") ~ "Tobolka",
      str_detect(drug_form, "injekční|Injekční|infuzní|Infuzní") ~ "Injekce / Infuze",
      str_detect(drug_form, "č[íi]pek|Čípek") ~ "Čípek",
      str_detect(drug_form, "perorální roztok|Perorální roztok") ~ "Perorální roztok",
      TRUE ~ "Other"
    )
  )

prescriptions <- prescriptions %>%
  mutate(
    year = year(presc_date)
  )

ggplot(prescriptions %>% filter(!is.na(presc_date), !is.na(drug_form)), aes(x = factor(year), fill = drug_form_std)) +
  geom_bar(position = "fill") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(
    title = "Poměr způsobu podání v čase",
    x = "Rok",
    y = "%",
    fill = "Způsob podání"
  ) +
  theme_minimal()

client_yearly_counts <- prescriptions %>%
  filter(!is.na(presc_date), !is.na(drug_form)) %>%
  group_by(client_id, year, drug_form_std) %>%
  summarise(presc_count = n(), .groups = "drop")

average_freq_per_year <- client_yearly_counts %>%
  group_by(year, drug_form_std) %>%
  summarise(avg_presc_per_client = mean(presc_count), .groups = "drop")

ggplot(average_freq_per_year, aes(x = factor(year), y = avg_presc_per_client, fill = drug_form_std)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(
    title = "Počet předpisů za rok na klienta dle způsobu podání",
    x = "Rok",
    y = "Průměrný počet předpisů na klienta",
    fill = "Způsob podání"
  ) +
  theme_minimal()

Výkyv v poměru předpisů ve druhé polovině roku 2022 je v případě této kohorty nejvyšší.

filtered_data <- prescriptions %>%
  mutate(
    drug_form_std = ifelse(is.na(drug_form_std), "Unknown", drug_form_std)
  ) %>%
  filter(
    drug_form_std %in% c("Injekce / Infuze", "Tableta"),
    !is.na(vaccinated),
    !is.na(presc_yearmonth)
  )

form_vax_counts <- filtered_data %>%
  group_by(presc_yearmonth, vaccinated, drug_form_std) %>%
  summarise(n = n(), .groups = "drop")

month_totals <- filtered_data %>%
  group_by(presc_yearmonth, vaccinated) %>%
  summarise(total = n(), .groups = "drop")

form_vax_ratios <- form_vax_counts %>%
  left_join(month_totals, by = c("presc_yearmonth", "vaccinated")) %>%
  mutate(ratio = n / total)

form_vax_ratios <- form_vax_ratios %>%
  mutate(group = paste(drug_form_std, ifelse(vaccinated == "Vakcinace", "Vakcinace", "Bez vakcinace"), sep = " - "))

ggplot(form_vax_ratios, aes(x = presc_yearmonth, y = ratio, color = group)) +
  geom_line(size = 0.6) +                   
  geom_smooth(se = FALSE, span = 0.3) +   
  scale_y_continuous(labels = percent_format(accuracy = 1)) +
  scale_x_date(date_labels = "%Y",         
               date_breaks = "1 year",
               expand = c(0.01, 0)) +
  labs(
    title = "Poměr počtu předpisů podle způsobu podání a vakcinace",
    x = "Datum",
    y = "Poměr počtu předpisů",
    color = "Forma podání + Vakcinace"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0, hjust = 0.5))

8.2 Krátkodobý efekt

8.2.1 Prvopředpisy

8.2.1.0.1 SCCS vakcinovaných

Riziko prvopředpisů u očkovaných je v zásadě nižší oproti baseline. Statisticky významné výsledky jsou však jen v časovém období 30-89 a 180-269.

combined_vaccinated_prescribed <- dbGetQuery(con,
                                    "
                                    SELECT
                                    v.client_id,
                                    v.date AS vaccination_date,
                                    p.date AS event_date,
                                    'CPZP' as ins_company
                                    FROM cpzp_vaccinations v
                                    LEFT JOIN cpzp_clients c ON v.client_id = c.client_id
                                    LEFT JOIN cpzp_ins_start_end s ON c.client_id = s.client_id
                                    INNER JOIN cpzp_prescriptions p ON c.client_id = p.client_id
                                                                    AND p.presc_order = 1
                                                                    AND p.date BETWEEN v.date - INTERVAL '1 year' AND v.date + INTERVAL '1 year'
                                    WHERE v.event_order = 1
                                    AND s.ins_start_date < v.date - INTERVAL '1 year'
                                    AND s.ins_end_date > v.date + INTERVAL '1 year'
                                    AND v.date BETWEEN DATE '2021-05-21' - INTERVAL '7 days'
                                                       AND DATE '2021-05-21' + INTERVAL '7 days'
                                    AND 2021 - c.birthyear BETWEEN 40 AND 59
                                    
                                    UNION
                                    
                                    SELECT
                                    v.client_id,
                                    v.date AS vaccination_date,
                                    p.date AS event_date,
                                    'OZP' AS ins_company
                                    FROM OZP_vaccinations v
                                    LEFT JOIN ozp_clients c ON v.client_id = c.client_id
                                    LEFT JOIN ozp_ins_start_end s ON c.client_id = s.client_id
                                    INNER JOIN ozp_prescriptions p ON c.client_id = p.client_id
                                                                    AND p.presc_order = 1
                                                                    AND p.date BETWEEN v.date - INTERVAL '1 year' AND v.date + INTERVAL '1 year'
                                    WHERE v.event_order = 1
                                    AND s.ins_start_date < v.date - INTERVAL '1 year'
                                    AND s.ins_end_date > v.date + INTERVAL '1 year'
                                    AND v.date BETWEEN DATE '2021-05-21' - INTERVAL '7 days'
                                                       AND DATE '2021-05-21' + INTERVAL '7 days'
                                    AND 2021 - c.birthyear BETWEEN 40 AND 59
                                    ") %>%
  mutate(client_id = row_number()) %>%
  transmute(
    person_id = client_id,
    exposure_start_date = vaccination_date,
    event_date,
    observation_start_date = vaccination_date - days(365),
    observation_end_date = vaccination_date + days(365)
  )

sccs_vaccinated <- combined_vaccinated_prescribed %>%
  mutate(
    event_day = as.integer(event_date - observation_start_date),
    exposure_start_day = as.integer(exposure_start_date - observation_start_date),
    observation_start_day = 0,
    observation_end_day = as.integer(observation_end_date - observation_start_date)
  )

sccs_result_vaccinated <- standardsccs(
  formula = event_day ~ exposure_start_day,
  indiv = person_id,
  astart = observation_start_day,
  aend = observation_end_day,
  aevent = event_day,
  adrug = exposure_start_day,
  aedrug = exposure_start_day + 365,
  expogrp = list(c(0, 30, 90, 180, 270)),
  data = sccs_vaccinated
)

cat("```\n", paste(capture.output(print(sccs_result_vaccinated)), collapse = "\n"), "\n```")
 Call:
coxph(formula = Surv(rep(1, 34968L), event) ~ exposure_start_day + 
    strata(indivL) + offset(log(interval)), data = chopdat, method = "exact")

  n= 34968, number of events= 5828 

                        coef exp(coef) se(coef)      z Pr(>|z|)    
exposure_start_day1 -0.07643   0.92642  0.06796 -1.125  0.26074    
exposure_start_day2 -0.28572   0.75147  0.05452 -5.241  1.6e-07 ***
exposure_start_day3 -0.05380   0.94762  0.04154 -1.295  0.19525    
exposure_start_day4 -0.10989   0.89593  0.04250 -2.586  0.00972 ** 
exposure_start_day5 -0.04944   0.95176  0.04040 -1.224  0.22107    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                    exp(coef) exp(-coef) lower .95 upper .95
exposure_start_day1    0.9264      1.079    0.8109    1.0584
exposure_start_day2    0.7515      1.331    0.6753    0.8362
exposure_start_day3    0.9476      1.055    0.8735    1.0280
exposure_start_day4    0.8959      1.116    0.8243    0.9738
exposure_start_day5    0.9518      1.051    0.8793    1.0302

Concordance= 0.763  (se = 0.005 )
Likelihood ratio test= 32.93  on 5 df,   p=4e-06
Wald test            = 31.21  on 5 df,   p=9e-06
Score (logrank) test = 31.36  on 5 df,   p=8e-06
 
8.2.1.0.2 SCCS nevakcinovaných

Zde je zajímavá změna trendu, kdy riziko neočkovaných v období 180-269 mírně přesáhne baseline.

combined_unvaccinated_prescribed <- dbGetQuery(con,
                                               "
                                    SELECT
                                    p.client_id,
                                    DATE '2021-05-21' AS vaccination_date,
                                    p.date AS event_date,
                                    'CPZP' as ins_company
                                    FROM cpzp_prescriptions p
                                    LEFT JOIN cpzp_clients c ON P.client_id = c.client_id 
                                    LEFT JOIN cpzp_ins_start_end s ON c.client_id = s.client_id
                                    WHERE 1 = 1
                                    AND s.ins_start_date < DATE '2021-05-21'- INTERVAL '1 year'
                                    AND s.ins_end_date > DATE '2021-05-21' + INTERVAL '1 year'
                                    AND 2021 - c.birthyear BETWEEN 40 AND 59
                                    AND p.date BETWEEN DATE '2021-05-21' - INTERVAL '1 year' AND DATE '2021-05-21' + INTERVAL '1 year'
                                    AND p.presc_order = 1
                                    AND c.n_vacc = 0
                                    
                                    UNION
                                    
                                    SELECT
                                    p.client_id,
                                    DATE '2021-05-21' AS vaccination_date,
                                    p.date AS event_date,
                                    'OZP' as ins_company
                                    FROM ozp_prescriptions p
                                    LEFT JOIN ozp_clients c ON P.client_id = c.client_id 
                                    LEFT JOIN ozp_ins_start_end s ON c.client_id = s.client_id
                                    WHERE 1 = 1
                                    AND s.ins_start_date < DATE '2021-05-21' - INTERVAL '1 year'
                                    AND s.ins_end_date > DATE '2021-05-21' + INTERVAL '1 year'
                                    AND 2021 - c.birthyear BETWEEN 40 AND 59
                                    AND p.date BETWEEN DATE '2021-05-21' - INTERVAL '1 year' AND DATE '2021-05-21' + INTERVAL '1 year'
                                    AND p.presc_order = 1
                                    AND c.n_vacc = 0

                                               ") %>%
  mutate(client_id = row_number()) %>%
  transmute(
    person_id = client_id,
    exposure_start_date = vaccination_date,
    event_date,
    observation_start_date = vaccination_date - days(365),
    observation_end_date = vaccination_date + days(365)
  )

sccs_unvaccinated <- combined_unvaccinated_prescribed %>%
  mutate(
    event_day = as.integer(event_date - observation_start_date),
    exposure_start_day = as.integer(exposure_start_date - observation_start_date),
    observation_start_day = 0,
    observation_end_day = as.integer(observation_end_date - observation_start_date)
  )

sccs_result_unvaccinated <- standardsccs(
  formula = event_day ~ exposure_start_day,
  indiv = person_id,
  astart = observation_start_day,
  aend = observation_end_day,
  aevent = event_day,
  adrug = exposure_start_day,
  aedrug = exposure_start_day + 365,
  expogrp = list(c(0, 30, 90, 180, 270)),
  data = sccs_unvaccinated
)

cat("```\n", paste(capture.output(print(sccs_result_unvaccinated)), collapse = "\n"), "\n```")
 Call:
coxph(formula = Surv(rep(1, 45756L), event) ~ exposure_start_day + 
    strata(indivL) + offset(log(interval)), data = chopdat, method = "exact")

  n= 45756, number of events= 7626 

                        coef exp(coef) se(coef)      z Pr(>|z|)   
exposure_start_day1 -0.03966   0.96111  0.05977 -0.664  0.50696   
exposure_start_day2 -0.14062   0.86882  0.04575 -3.074  0.00211 **
exposure_start_day3  0.03271   1.03325  0.03589  0.911  0.36215   
exposure_start_day4  0.06801   1.07038  0.03539  1.922  0.05465 . 
exposure_start_day5 -0.03389   0.96668  0.03592 -0.943  0.34552   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                    exp(coef) exp(-coef) lower .95 upper .95
exposure_start_day1    0.9611     1.0405    0.8549    1.0806
exposure_start_day2    0.8688     1.1510    0.7943    0.9503
exposure_start_day3    1.0332     0.9678    0.9631    1.1086
exposure_start_day4    1.0704     0.9342    0.9986    1.1473
exposure_start_day5    0.9667     1.0345    0.9010    1.0372

Concordance= 0.749  (se = 0.004 )
Likelihood ratio test= 18.4  on 5 df,   p=0.002
Wald test            = 18.04  on 5 df,   p=0.003
Score (logrank) test = 18.07  on 5 df,   p=0.003
 

8.2.1.1 Porovnání rizik

Při porovnání rizik je všeobecně méně riziková skupina ta očkovaná (která ani v jednom případě nepřesáhne úroveň zvýšeného rizika oproti baseline). Statisticky významné rozdíly jsou v období 30-89 a 180-269, kdy jsou rizika očkovaných významně menší.

coef_unvax <- sccs_result_unvaccinated$coefficients
conf_unvax <- sccs_result_unvaccinated$conf.int

df_unvax <- data.frame(
  term = rownames(coef_unvax),
  logHR = coef_unvax[, "coef"],
  SE = coef_unvax[, "se(coef)"],
  HR = coef_unvax[, "exp(coef)"],
  conf.low = conf_unvax[, "lower .95"],
  conf.high = conf_unvax[, "upper .95"]
)

coef_vax <- sccs_result_vaccinated$coefficients
conf_vax <- sccs_result_vaccinated$conf.int

df_vax <- data.frame(
  term = rownames(coef_vax),
  logHR = coef_vax[, "coef"],
  SE = coef_vax[, "se(coef)"],
  HR = coef_vax[, "exp(coef)"],
  conf.low = conf_vax[, "lower .95"],
  conf.high = conf_vax[, "upper .95"]
)

df_unvax <- df_unvax %>% mutate(group = "Bez Vakcinace")
df_vax <- df_vax %>% mutate(group = "Vakcinace")

combined <- bind_rows(df_vax, df_unvax)

combined$period <- as.integer(gsub("exposure_start_day", "", combined$term))

start_days <- c(0, 30, 90, 180, 270)
end_days   <- c(29, 89, 179, 269, Inf)

period_labels <- ifelse(
  is.infinite(end_days),
  paste0(start_days, "+"),
  paste0(start_days, "-", end_days)
)

combined$period_label <- factor(combined$period,
                                levels = seq_along(period_labels),
                                labels = period_labels)

ggplot(combined, aes(x = period_label, y = HR, color = group)) +
  geom_point(position = position_dodge(width = 0.6), size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2,
                position = position_dodge(width = 0.6)) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  labs(
    x = "Období expozice (dny)",
    y = "Poměr rizik (HR)",
    title = "Porovnání SCCS: očkovaní vs neočkovaní"
  ) +
  theme_minimal()

vax_df <- combined %>%
  filter(group == "Vakcinace") %>%
  rename_with(~ paste0(., "_vax"), -c(period, period_label))

unvax_df <- combined %>%
  filter(group == "Bez Vakcinace") %>%
  rename_with(~ paste0(., "_unvax"), -c(period, period_label))

df_diff <- inner_join(vax_df, unvax_df, by = c("period", "period_label"))

df_diff <- df_diff %>%
  mutate(
    logHR_diff = logHR_vax - logHR_unvax,
    se_diff = sqrt(SE_vax^2 + SE_unvax^2),
    ci_low = logHR_diff - 1.96 * se_diff,
    ci_high = logHR_diff + 1.96 * se_diff
  )

ggplot(df_diff, aes(x = period_label, y = logHR_diff)) +
  geom_point(color = "blue", size = 3) +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.2, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "Rozdíl v log(HR) (Očkovaní - Neočkovaní)",
    x = "Období expozice",
    y = "Rozdíl v log(HR)"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

common_terms <- intersect(df_vax$term, df_unvax$term)

did_results <- data.frame(
  term = character(),
  period_label = character(),
  logHR_vax = numeric(),
  logHR_unvax = numeric(),
  DiD_logHR = numeric(),
  DiD_SE = numeric(),
  Z = numeric(),
  P_value = numeric(),
  stringsAsFactors = FALSE
)

for (term in common_terms) {
  logHR_vax <- df_vax %>% filter(term == !!term) %>% pull(logHR)
  SE_vax    <- df_vax %>% filter(term == !!term) %>% pull(SE)
  
  logHR_unvax <- df_unvax %>% filter(term == !!term) %>% pull(logHR)
  SE_unvax    <- df_unvax %>% filter(term == !!term) %>% pull(SE)
  
  DiD_logHR <- logHR_vax - logHR_unvax
  DiD_SE <- sqrt(SE_vax^2 + SE_unvax^2)
  Z <- DiD_logHR / DiD_SE
  P_value <- 2 * pnorm(-abs(Z))
  
  period_idx <- as.integer(gsub("exposure_start_day", "", term))
  period_label <- period_labels[period_idx]
  
  did_results <- rbind(did_results, data.frame(
    term = term,
    period_label = period_label,
    logHR_vax = logHR_vax,
    logHR_unvax = logHR_unvax,
    DiD_logHR = DiD_logHR,
    DiD_SE = DiD_SE,
    Z = Z,
    P_value = P_value
  ))
}

did_results %>%
  transmute(Období = period_label, `log(HR) Vakcinace` = logHR_vax, `log(HR) Bez vakcinace` = logHR_unvax, `DiD log(HR)` = DiD_logHR, p_value = P_value) %>%
  kable(format = "html", caption = "Porovnání rizik očkovaných a neočkovaných (SCCS)") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání rizik očkovaných a neočkovaných (SCCS)
Období log(HR) Vakcinace log(HR) Bez vakcinace DiD log(HR) p_value
0-29 -0.0764318 -0.0396637 -0.0367680 0.6845618
30-89 -0.2857232 -0.1406190 -0.1451042 0.0414702
90-179 -0.0538000 0.0327083 -0.0865082 0.1150647
180-269 -0.1098894 0.0680120 -0.1779014 0.0012968
270+ -0.0494389 -0.0338858 -0.0155531 0.7735826

8.2.2 Předpisy celkem

8.2.2.1 SCCS vakcinovaných

Díky opravdu velkému množství dat jsou veškeré výsledky statisticky významné. Riziko jako takové kolísá nejprve snížením pod úroveň baseline a poté postupným navyšováním s vrcholem na konci sledovaného období.

combined_vaccinated_prescribed <- dbGetQuery(con,
                                    "
                                    SELECT
                                    v.client_id,
                                    v.date AS vaccination_date,
                                    p.date AS event_date,
                                    d.drug_form,
                                    'CPZP' as ins_company
                                    FROM cpzp_vaccinations v
                                    LEFT JOIN cpzp_clients c ON v.client_id = c.client_id
                                    LEFT JOIN cpzp_ins_start_end s ON c.client_id = s.client_id
                                    INNER JOIN cpzp_prescriptions p ON c.client_id = p.client_id
                                                                    AND p.date BETWEEN v.date - INTERVAL '1 year' AND v.date + INTERVAL '1 year'
                                    LEFT JOIN cpzp_drugs_catalog d ON p.presc_id  = d.presc_id 
                                    WHERE v.event_order = 1
                                    AND s.ins_start_date < v.date - INTERVAL '1 year'
                                    AND s.ins_end_date > v.date + INTERVAL '1 year'
                                    AND v.date BETWEEN DATE '2021-05-21' - INTERVAL '7 days'
                                                       AND DATE '2021-05-21' + INTERVAL '7 days'
                                    AND 2021 - c.birthyear BETWEEN 40 AND 59
                                    
                                    UNION
                                    
                                    SELECT
                                    v.client_id,
                                    v.date AS vaccination_date,
                                    p.date AS event_date,
                                    d.drug_form,
                                    'OZP' AS ins_company
                                    FROM OZP_vaccinations v
                                    LEFT JOIN ozp_clients c ON v.client_id = c.client_id
                                    LEFT JOIN ozp_ins_start_end s ON c.client_id = s.client_id
                                    INNER JOIN ozp_prescriptions p ON c.client_id = p.client_id
                                                                    AND p.date BETWEEN v.date - INTERVAL '1 year' AND v.date + INTERVAL '1 year'
                                    LEFT JOIN ozp_drugs_catalog d ON p.presc_id  = d.presc_id 
                                    WHERE v.event_order = 1
                                    AND s.ins_start_date < v.date - INTERVAL '1 year'
                                    AND s.ins_end_date > v.date + INTERVAL '1 year'
                                    AND v.date BETWEEN DATE '2021-05-21' - INTERVAL '7 days'
                                                       AND DATE '2021-05-21' + INTERVAL '7 days'
                                    AND 2021 - c.birthyear BETWEEN 40 AND 59
                                    ") %>%
  mutate(client_id = as.numeric(factor(client_id)),
         drug_form = tolower(drug_form),
         inj_or_inf = str_detect(drug_form, "injekční|infuzní"))

injectable_infusional <- combined_vaccinated_prescribed %>%
  filter(inj_or_inf) %>%
  arrange(client_id, event_date) %>%
  group_by(client_id) %>%
  mutate(
    day_diff = as.numeric(difftime(event_date, lag(event_date, default = first(event_date)), units = "days")),
    new_window = is.na(day_diff) | day_diff > 10,
    window_id = cumsum(new_window)
  ) %>%
  group_by(client_id, window_id) %>%
  slice(1) %>%
  ungroup()

other_drugs <- combined_vaccinated_prescribed %>%
  filter(!inj_or_inf)

combined_vaccinated_prescribed <- bind_rows(injectable_infusional, other_drugs) %>%
  arrange(client_id, event_date) %>%
  transmute(
    person_id = client_id,
    exposure_start_date = vaccination_date,
    event_date,
    observation_start_date = vaccination_date - days(365),
    observation_end_date = vaccination_date + days(365)
  )

sccs_vaccinated <- combined_vaccinated_prescribed %>%
  mutate(
    event_day = as.integer(event_date - observation_start_date),
    exposure_start_day = as.integer(exposure_start_date - observation_start_date),
    observation_start_day = 0,
    observation_end_day = as.integer(observation_end_date - observation_start_date)
  )

sccs_result_vaccinated <- standardsccs(
  formula = event_day ~ exposure_start_day,
  indiv = person_id,
  astart = observation_start_day,
  aend = observation_end_day,
  aevent = event_day,
  adrug = exposure_start_day,
  aedrug = exposure_start_day + 365,
  expogrp = list(c(0, 30, 90, 180, 270)),
  data = sccs_vaccinated
)

cat("```\n", paste(capture.output(print(sccs_result_vaccinated)), collapse = "\n"), "\n```")
 Call:
coxph(formula = Surv(rep(1, 183642L), event) ~ exposure_start_day + 
    strata(indivL) + offset(log(interval)), data = chopdat, method = "exact")

  n= 183642, number of events= 30607 

                        coef exp(coef) se(coef)      z Pr(>|z|)    
exposure_start_day1  0.07770   1.08080  0.02882  2.696 0.007011 ** 
exposure_start_day2 -0.14287   0.86687  0.02331 -6.130 8.82e-10 ***
exposure_start_day3  0.06440   1.06651  0.01804  3.570 0.000357 ***
exposure_start_day4  0.08328   1.08685  0.01791  4.651 3.30e-06 ***
exposure_start_day5  0.16960   1.18483  0.01689 10.043  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                    exp(coef) exp(-coef) lower .95 upper .95
exposure_start_day1    1.0808     0.9252    1.0214    1.1436
exposure_start_day2    0.8669     1.1536    0.8282    0.9074
exposure_start_day3    1.0665     0.9376    1.0295    1.1049
exposure_start_day4    1.0868     0.9201    1.0494    1.1257
exposure_start_day5    1.1848     0.8440    1.1463    1.2247

Concordance= 0.744  (se = 0.002 )
Likelihood ratio test= 183.4  on 5 df,   p=<2e-16
Wald test            = 182.6  on 5 df,   p=<2e-16
Score (logrank) test = 183.2  on 5 df,   p=<2e-16
 

8.2.2.2 SCCS nevakcinovaných

Trend neočkovaných v zásadě kopíruje trend vakcinované skupiny a to i v ohledu významnosti výsledků.

combined_unvaccinated_prescribed <- dbGetQuery(con,
                                               "
                                    SELECT
                                    p.client_id,
                                    DATE '2021-05-21' AS vaccination_date,
                                    p.date AS event_date,
                                    d.drug_form,
                                    'CPZP' as ins_company
                                    FROM cpzp_prescriptions p
                                    LEFT JOIN cpzp_clients c ON P.client_id = c.client_id 
                                    LEFT JOIN cpzp_ins_start_end s ON c.client_id = s.client_id
                                    LEFT JOIN cpzp_drugs_catalog d ON p.presc_id  = d.presc_id 
                                    WHERE 1 = 1
                                    AND s.ins_start_date < DATE '2021-05-21' - INTERVAL '1 year'
                                    AND s.ins_end_date > DATE '2021-05-21' + INTERVAL '1 year'
                                    AND 2021 - c.birthyear BETWEEN 40 AND 59
                                    AND p.date BETWEEN DATE '2021-05-21' - INTERVAL '1 year' AND DATE '2021-05-21' + INTERVAL '1 year'
                                    AND c.n_vacc = 0
                                    
                                    UNION
                                    
                                    SELECT
                                    p.client_id,
                                    DATE '2021-05-21' AS vaccination_date,
                                    p.date AS event_date,
                                    d.drug_form,
                                    'OZP' as ins_company
                                    FROM ozp_prescriptions p
                                    LEFT JOIN ozp_clients c ON P.client_id = c.client_id 
                                    LEFT JOIN ozp_ins_start_end s ON c.client_id = s.client_id
                                    LEFT JOIN ozp_drugs_catalog d ON p.presc_id  = d.presc_id 
                                    WHERE 1 = 1
                                    AND s.ins_start_date < DATE '2021-05-21' - INTERVAL '1 year'
                                    AND s.ins_end_date > DATE '2021-05-21' + INTERVAL '1 year'
                                    AND 2021 - c.birthyear BETWEEN 40 AND 59
                                    AND p.date BETWEEN DATE '2021-05-21' - INTERVAL '1 year' AND DATE '2021-05-21' + INTERVAL '1 year'
                                    AND c.n_vacc = 0

                                               ") %>%
  mutate(client_id = as.numeric(factor(client_id)),
         drug_form = tolower(drug_form),
         inj_or_inf = str_detect(drug_form, "injekční|infuzní"))

injectable_infusional <- combined_unvaccinated_prescribed %>%
  filter(inj_or_inf) %>%
  arrange(client_id, event_date) %>%
  group_by(client_id) %>%
  mutate(
    day_diff = as.numeric(difftime(event_date, lag(event_date, default = first(event_date)), units = "days")),
    new_window = is.na(day_diff) | day_diff > 10,
    window_id = cumsum(new_window)
  ) %>%
  group_by(client_id, window_id) %>%
  slice(1) %>%
  ungroup()

other_drugs <- combined_unvaccinated_prescribed %>%
  filter(!inj_or_inf)

combined_unvaccinated_prescribed <- bind_rows(injectable_infusional, other_drugs) %>%
  arrange(client_id, event_date) %>%
  transmute(
    person_id = client_id,
    exposure_start_date = vaccination_date,
    event_date,
    observation_start_date = vaccination_date - days(365),
    observation_end_date = vaccination_date + days(365)
  )

sccs_unvaccinated <- combined_unvaccinated_prescribed %>%
  mutate(
    event_day = as.integer(event_date - observation_start_date),
    exposure_start_day = as.integer(exposure_start_date - observation_start_date),
    observation_start_day = 0,
    observation_end_day = as.integer(observation_end_date - observation_start_date)
  )

sccs_result_unvaccinated <- standardsccs(
  formula = event_day ~ exposure_start_day,
  indiv = person_id,
  astart = observation_start_day,
  aend = observation_end_day,
  aevent = event_day,
  adrug = exposure_start_day,
  aedrug = exposure_start_day + 365,
  expogrp = list(c(0, 30, 90, 180, 270)),
  data = sccs_unvaccinated
)

cat("```\n", paste(capture.output(print(sccs_result_unvaccinated)), collapse = "\n"), "\n```")
 Call:
coxph(formula = Surv(rep(1, 236910L), event) ~ exposure_start_day + 
    strata(indivL) + offset(log(interval)), data = chopdat, method = "exact")

  n= 236910, number of events= 39485 

                        coef exp(coef) se(coef)      z Pr(>|z|)    
exposure_start_day1  0.04053   1.04136  0.02593  1.563   0.1181    
exposure_start_day2 -0.06283   0.93910  0.01991 -3.155   0.0016 ** 
exposure_start_day3  0.06723   1.06954  0.01594  4.217 2.47e-05 ***
exposure_start_day4  0.11704   1.12417  0.01563  7.487 7.02e-14 ***
exposure_start_day5  0.17413   1.19021  0.01491 11.675  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                    exp(coef) exp(-coef) lower .95 upper .95
exposure_start_day1    1.0414     0.9603    0.9898    1.0957
exposure_start_day2    0.9391     1.0649    0.9032    0.9765
exposure_start_day3    1.0695     0.9350    1.0366    1.1035
exposure_start_day4    1.1242     0.8895    1.0902    1.1591
exposure_start_day5    1.1902     0.8402    1.1559    1.2255

Concordance= 0.742  (se = 0.002 )
Likelihood ratio test= 197.1  on 5 df,   p=<2e-16
Wald test            = 199.9  on 5 df,   p=<2e-16
Score (logrank) test = 200.3  on 5 df,   p=<2e-16
 

8.2.2.3 Porovnání rizik

Vzhledem k tomu, že trendy jsou v zásadě srovnatelné, není možné po většinu sledovaného období jednoznačně vyvodit vyšší riziko pro jednu nebo druhou skupinu. Jediný rozdíl je v období 30-89, kdy je riziko významně nižší pro očkovanou populaci.

coef_unvax <- sccs_result_unvaccinated$coefficients
conf_unvax <- sccs_result_unvaccinated$conf.int

df_unvax <- data.frame(
  term = rownames(coef_unvax),
  logHR = coef_unvax[, "coef"],
  SE = coef_unvax[, "se(coef)"],
  HR = coef_unvax[, "exp(coef)"],
  conf.low = conf_unvax[, "lower .95"],
  conf.high = conf_unvax[, "upper .95"]
)

coef_vax <- sccs_result_vaccinated$coefficients
conf_vax <- sccs_result_vaccinated$conf.int

df_vax <- data.frame(
  term = rownames(coef_vax),
  logHR = coef_vax[, "coef"],
  SE = coef_vax[, "se(coef)"],
  HR = coef_vax[, "exp(coef)"],
  conf.low = conf_vax[, "lower .95"],
  conf.high = conf_vax[, "upper .95"]
)

df_unvax <- df_unvax %>% mutate(group = "Unvaccinated")
df_vax <- df_vax %>% mutate(group = "Vaccinated")

combined <- bind_rows(df_vax, df_unvax)

combined$period <- as.integer(gsub("exposure_start_day", "", combined$term))

start_days <- c(0, 30, 90, 180, 270)
end_days   <- c(29, 89, 179, 269, Inf)

period_labels <- ifelse(
  is.infinite(end_days),
  paste0(start_days, "+"),
  paste0(start_days, "-", end_days)
)

combined$period_label <- factor(combined$period,
                                levels = seq_along(period_labels),
                                labels = period_labels)

ggplot(combined, aes(x = period_label, y = HR, color = group)) +
  geom_point(position = position_dodge(width = 0.6), size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2,
                position = position_dodge(width = 0.6)) +
  geom_hline(yintercept = 1, linetype = "dashed") +
  labs(
    x = "Období expozice (dny)",
    y = "Poměr rizik (HR)",
    title = "Porovnání SCCS: očkovaní vs neočkovaní"
  ) +
  theme_minimal()

vax_df <- combined %>%
  filter(group == "Vaccinated") %>%
  rename_with(~ paste0(., "_vax"), -c(period, period_label))

unvax_df <- combined %>%
  filter(group == "Unvaccinated") %>%
  rename_with(~ paste0(., "_unvax"), -c(period, period_label))

df_diff <- inner_join(vax_df, unvax_df, by = c("period", "period_label"))

df_diff <- df_diff %>%
  mutate(
    logHR_diff = logHR_vax - logHR_unvax,
    se_diff = sqrt(SE_vax^2 + SE_unvax^2),
    ci_low = logHR_diff - 1.96 * se_diff,
    ci_high = logHR_diff + 1.96 * se_diff
  )

ggplot(df_diff, aes(x = period_label, y = logHR_diff)) +
  geom_point(color = "blue", size = 3) +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = 0.2, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "Rozdíl v log(HR) (Očkovaní - Neočkovaní)",
    x = "Období expozice",
    y = "Rozdíl v log(HR)"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

common_terms <- intersect(df_vax$term, df_unvax$term)

did_results <- data.frame(
  term = character(),
  period_label = character(),
  logHR_vax = numeric(),
  logHR_unvax = numeric(),
  DiD_logHR = numeric(),
  DiD_SE = numeric(),
  Z = numeric(),
  P_value = numeric(),
  stringsAsFactors = FALSE
)

for (term in common_terms) {
  logHR_vax <- df_vax %>% filter(term == !!term) %>% pull(logHR)
  SE_vax    <- df_vax %>% filter(term == !!term) %>% pull(SE)
  
  logHR_unvax <- df_unvax %>% filter(term == !!term) %>% pull(logHR)
  SE_unvax    <- df_unvax %>% filter(term == !!term) %>% pull(SE)
  
  DiD_logHR <- logHR_vax - logHR_unvax
  DiD_SE <- sqrt(SE_vax^2 + SE_unvax^2)
  Z <- DiD_logHR / DiD_SE
  P_value <- 2 * pnorm(-abs(Z))
  
  period_idx <- as.integer(gsub("exposure_start_day", "", term))
  period_label <- period_labels[period_idx]
  
  did_results <- rbind(did_results, data.frame(
    term = term,
    period_label = period_label,
    logHR_vax = logHR_vax,
    logHR_unvax = logHR_unvax,
    DiD_logHR = DiD_logHR,
    DiD_SE = DiD_SE,
    Z = Z,
    P_value = P_value
  ))
}

did_results %>%
  transmute(Období = period_label, `log(HR) Vakcinace` = logHR_vax, `log(HR) Bez vakcinace` = logHR_unvax, `DiD log(HR)` = DiD_logHR, p_value = P_value) %>%
  kable(format = "html", caption = "Porovnání rizik (SCCS) očkovaných a neočkovaných") %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed"),
    full_width = TRUE
  )
Porovnání rizik (SCCS) očkovaných a neočkovaných
Období log(HR) Vakcinace log(HR) Bez vakcinace DiD log(HR) p_value
0-29 0.0777049 0.0405313 0.0371735 0.3376538
30-89 -0.1428676 -0.0628347 -0.0800329 0.0090377
90-179 0.0643952 0.0672268 -0.0028315 0.9063690
180-269 0.0832830 0.1170410 -0.0337580 0.1555280
270+ 0.1696027 0.1741257 -0.0045229 0.8408985

8.3 Shrnutí

Oproti ostatním kohortám dochází ve věkové kategorii 40-59 k mírnému propadu v prvopředpisovosti, kde mezi jednotlivými skupinami není významný rozdíl. Naopak měsíční multipředpisovost vykazuje nepatrný nárůst se statisticky signifikantním rozdílem mezi skupinami (konkrétně téměř nulový nárůst u očkovaných a zhruba 7% nárůst u neočkované populace.).

Trendy ve formě podání jsou téměř konstantní s jedinou výjimkou opět ve druhé polovině roku 2021.

Z hlediska krátkodobých rizik se trendy obou skupin víceméně kopírují s několika málo významnými odchylkami, kde je nižší riziko vypočtené u očkovaných.

9 Závěr

Všeobecně je možné konstatovat, že riziko předpisu kortikoidsteroidů během covidu a v období po covidu rostlo a neustále roste. Co se týče vlivu očkování na tento trend, tak ve většině případů není možné zhodnotit rozdíly mezi skupinami jako významné. V několika případech, kdy se riziko liší, je vakcinovaná populace vyhodnocená jako méně riziková s jedinou výjimkou u rizik celkového počtu předpisů 270-365 dní po vakcinaci pro kohortu ve věku 16-39 let.

V zásadě je také možné říct, že rizika se týkala zejména mladší populace. V případě starší populace sledujeme spíše konstantní rizika a nebo v některých případech dokonce nepatrné snížení.

Ve zkratce tedy:

Změnila se v období pandemie nějak celková spotřeba kortikoidů? Ano. Nejprve došlo k očekávatelnému propadu a poté došlo k navýšení, které přesahuje trend nastavený předcovidovým obdobím.

Jsou trendy týkající se spotřeby kortikoidů stejné u očkované i neočkované populace? Chránily tedy vakcíny před nárůstem spotřeby kortikoidů, neměly žádný efekt, nebo dokonce škodily? Trendy jsou v zásadě srovnatelné a v případech, kdy tomu tak není, je vakcinovaná populace vyhodnocená jako s nižším rizikem.

Je rozdíl mezi jednotlivými věkovými skupinami? Ano. K největším rozdílům oproti celkovému trendu dochází u mladší populace. Je však otázka zdali to není způsobené prostým stárnutím dané kohorty.

Je důležité zmínit, že samotné zvýšení počtu předpisů nemusí nutně znamenat zhoršení situace. Vzhledem k tomu, že se kortikoidy předepisují na celou škálu různých problémů a každý problém má rozdílné dávkovací schéma, může jednoduše docházet i k nárůstu diagnóz s nutností častějšího (nebo sezónního) dávkování. Vzhledem k naší odbornosti by ale bylo případné zhodnocení spíše spekulací.

Zajímavé bylo pozorovat výrazný rozdíl v poměru předepsaných injekcí / infuzí na tablety u neočkované populace ve druhé polovině roku 2022. Došlo zde totiž k výraznému snížení počtu injekcí / infuzí a zároveň i zvýšení tablet. Zdali se jedná o výsledky covidových nařízení nebo jiné důvody opět přenecháme erudovanějším osobám.

Během zkoumání bylo (přirozeně) nalezeno mnoho slepých uliček, ale i zajímavých dat, která nejsou v tomto textu z kapacitních důvodů uvedena. Zejména data týkající se multipředpisovosti mohou ukázat mnoho dalších směrů k bádání. Pro příklad: průměrné intervaly mezi jednotlivými předpisy na omezeném intervalu 0-365, míra rizika druhého předpisu po prvním předpisu v čase, vliv incidence COVID-19 na předpisovost kortikoidů, průměrný počet předpisů na multipředepisovaného klienta, změny v sezónnosti předpisů a mnoho dalších.